Simulation and inference for a stochastic process with a random time of event in rstan

I would like to model a continuous-time (partly) stochastic process X(t) that has a jump at a random time point \tau and that is not directly observed.

The process is deterministic up to \tau

X(t) \equiv 0 \qquad \text{ for } t < \tau,

at \tau, it jumps to a random value m_0

X(\tau) = m_0,

and after \tau, it evolves as a Markov process with a transition probability that depends on a parameter \theta:

\left\lbrace X(t)\, {\large\vert}\, X(s) = x_s\right\rbrace \sim \mathcal{N}_{>0} \left(x_s - \theta \cdot x_s\cdot (t-s), \,\sqrt{\theta \cdot x_s\cdot (t-s)}\right)\qquad \text{ for } \tau\leq s < t.

I assume multiplicative measurement errors for the observations y_1,\dots,y_M at time points t_1,\dots, t_M:

\log(y_k)\sim\mathcal{N}(\log(X(t_k) + 1),\, \sigma)\qquad \text{ for } k=1,\dots M.

Moreover, I assume the following priors

\tau \sim \mathcal{Unif}(0,10), \quad m_0 \sim \mathcal{N}_{>0}(300,300), \quad \theta \sim \mathcal{N}_{>0}(0,5), \quad \sigma\sim\mathcal{Unif}(0.001,2),

where \mathcal{N}_{>0} denotes a truncated normal distribution. Note that the random time point \tau does not need to be one of the observed time point.

I am able to simulate the model in rstan with simulate_stoch_exp_decay_with_random_release_time_fix_param.stan (2.1 KB) :

stansim <- stan(file = "simulate_stoch_exp_decay_with_random_release_time_fix_param.stan", 
                data = list(M=181, time_points=seq(from=0, to=30, length.out=181), 
                            tau=1.1, m0=200, theta=.11, sigma=.01), 
                algorithm = "Fixed_param", seed = 4242, iter = 1, chains = 1)
y_obs <- extract(stansim, pars = "y_obs")[[1]][1,]
write.table(y_obs, file = "y_obs.txt", col.names=FALSE, row.names=FALSE)

Now, I would like to infer the four parameters \tau, m_0, \theta, \sigma from the simulated dataset y_obs.txt (3.0 KB).

Among others, I’ve tried the following Stan Code (stanmodel_stoch_exp_decay_with_random_release_time_v1.stan (2.1 KB) ):

data {
  int<lower=0> M;
  vector<lower=0>[M] y_obs; // observations
  vector<lower=0>[M] time_points; // time points of the observations

// we assume log-normally distributed noise and therefore log-transform the data
transformed data{
  vector[M] log_y_obs;
  log_y_obs = log(y_obs);

parameters {
  real<lower=0, upper=10> tau; // time point at which initial amount m0 is released
  real<lower=0> m0; // initial amount released at time tau
  real<lower=0> theta; // decay rate
  real<lower=0.001, upper=2> sigma; // measurement noise parameter
  real<lower=0> x_k[M];// random process states

transformed parameters{
  real<lower=0> x[M]; // process states at observed time points
  for (t in 1:M) {
    if (time_points[t] < tau ){
      x[t] = 0; // process states before tau are deterministic
      if(time_points[t] == tau ){// tau is an observed time point
        x[t] = m0;
      }else{ // time_points[t] > tau
      x[t] = x_k[t];

model {
  // log-likelihood part 1: prod. of transition probabilities of the diffusion process
  int first; // variable to check whether the first transition after tau is yet to come
  first = 1;
  for (t in 1:M) {
    if(time_points[t] >= tau ){
      if(time_points[t] == tau ){ // tau is an observed time point
        first = 0;
      }else{// time_points[t] > tau
        if(first){// transtition prob. between tau and the first observed time point after tau
          x[t] ~ normal(m0 - theta * m0 * (time_points[t] - tau),
                      sqrt(theta * m0 * (time_points[t] - tau)));
          first = 0;
        }else{// transitions prob. between observed time points after tau
          x[t] ~ normal(x[t-1] - theta * x[t-1] * (time_points[t] - time_points[t-1]),
                        sqrt(theta * x[t-1] * (time_points[t] - time_points[t-1])));
  // log-likelihood part 2: multiplicative normal noise
  for (i in 1:M){
    log_y_obs[i] ~ normal(log(x[i] + 1), sigma);
  // priors
  theta ~ normal(0, 5);
  m0 ~ normal(300,300);

and in R:

y_obs <- read.table("y_obs.txt")[[1]]
stanfit1 <- stan(file = "stanmodel_stoch_exp_decay_with_random_release_time_v1.stan", 
                data = list(M=181, time_points=seq(from=0, to=30, length.out=181), 
                seed = 42424)

This results in mostly divergent transitions (increasing adapt_delta did not help) and the chains just get stuck:

I’m suspecting that sampling with Stan is simply not possible for this kind of model due to the fact that the components of X would have to switch between a deterministic and a stochastic state depending on the value of \tau. But I am not completely sure. Maybe there is some technique that I have not come across yet.

Does anyone have an idea on how to resolves this issue or can confirm my conjecture that sampling with Stan is not possible?