Ito process as numerical solution of stochastic differential equation

Using fitting trivial model dX= \theta dW as example, what I’d like to do in Stan is like this

data {
  int N;                        /* nb. of data points */
  int M;                        /* nb. of auxiliary data between each pair of data points */
  real<lower=0> T;              /* end time */
  vector[N] scal_w;             /* observed data X */
}

transformed data {
  real dt = T/(N*M);            /* fixed time step */
}
parameters {
  real<lower=0.0> sigma;        /* data simulated from given sigma */
  real<lower=0.0> eps;          /* obs error */
  real sn[M*N];                 /* Wiener process' std normal steps */
}

model {
  vector[M] x;
  real x0;
  
  sigma ~ normal(1.0, 0.5);
  eps ~ cauchy(0, 1);
  sn ~ normal(0, 1);

  x0 = 0.0;
  for (i in 1:N) {
                                /* nested loop equivalent to SDE integrator */
                                /* that outputs x, given sn, dt, and sigma */
    for (j in 1:M) {
      x[j] = (j == 1 ? x0 : x[j - 1]) + sigma * sqrt(dt) * sn[(i - 1) * M + j];
    }
    scal_w[i] ~ normal(x[M], eps);
    x0 = x[M];
  }
}

as proposed in http://apps.olin.wustl.edu/faculty/chib/papers/elerianchibshephard2001.pdf. The idea is to generate missing/auxiliary data between observations, with SDE integrator describing the recursive conditional likelihood p(X_{i+1}|X_i, W_{i+1}, \theta). This is a full-posterior approach compared with filter-based approaches.