The Merton Jump-Diffusion Model
Merton (1976) [1] jump-diffusion model is an extension to the Geometric Brownian Motion model, with the underlying asset exhibit jumps in addition to have continuous diffusion paths. The asset price evolves as 
dS t =μS t dt + σS t  dW t  + (η−1)dq  
where μ is the drift rate, σ the volatility of S  and dW t   the increment of a Wiener process. The independent Poisson process dq  has a value of 1  with probability λdt  and 0 otherwise. 

The following is a ThetaML implementation of the Merton Jump-Diffusion Model. 
Merton jump diffusion Process simulation
Implementation as ThetaML
model JumpDiffusionProcess
%Jump diffusion model for stock price movement
%Example input parameters:
%S0 = 100, r = 0.03, sigma = 0.2, lambda = 1
%nu = -0.1, delta = 0.05, t = 1, ht = 1/12
    import S0      "Initial stock price"
    import r       "Risk-free interest rate"
    import sigma   "Stock price volatility"
    import lambda  "Jump intensity"
    import nu      "Mean of jumps"
    import delta   "Std. of jumps"
    import t       "Simulation time horizon"
    import ht      "Discretization time step"
    export S       "Simulated stock prices"
 
    S = S0         %initialize stock prices
 
    loop t/ht      %loop t/ht times
        theta ht   %theta advances time by ht units
        P = poissrnd(lambda*ht,size(rand()))
        U = exp(P*nu + sqrt(P)*delta*randn())
        S = S * exp((r-lambda*(exp(nu+0.5*delta^2)-1)-0.5*sigma^2)*ht +      
              sigma*sqrt(ht)*randn())*U
    end
end

function state = Theta(dt, state)
 
  if nargin == 0
    state.S.value = 100;
    state.EUR.value = 1;
  else
    r = 0.05;
 
    mu1 = -0.2;
    lambda = 0.1;
    sigma = 0.2;
    gamma = 0.45;
 
    mu = r -(exp(mu1 + 0.5*gamma^2) - 1) * lambda;
 
    max_dt = lambda;
    n = max(dt/max_dt, 1);
 
    sub_dt = dt/n;
    dq = zeros(size(state.S));
    for i = 1 : n
      jumpsize = mu1 + gamma*randn(size(state.S));%
      u = rand(size(state.S));
      dq(u<lambda*sub_dt) = dq(u<lambda*sub_dt) + jumpsize(u<lambda*sub_dt);
    end
 
    state.S = state.S .* exp((mu-0.5*sigma^2)*dt + sqrt(dt)*sigma*randn(size(state.S))+dq);
    state.EUR = state.EUR* exp(-r*dt);
  end
end