Iâ€™m assuming you need y(t-\tau), not just (t-\tau), as in a time-dependent forcing term, e.g. \frac{dy}{dt}=f(y(t-\tau)) = sin(t-\tau) g(y), since that is normally feasible with ODE solvers.

I think the solution above is the best if your model can be discretized and solved with these simple differences instead of a more sophisticated solver (I often have to implement discrete stochastic versions of ODE models and compare their output to the mean field approximations and I rarely have problems, but of course that is model and parameterization-dependent). If that holds it would be actually what you are asking for, just being a â€śdelay difference equationâ€ť instead.

As an alternative, you may be able to formulate your model with a series of compartments before the component y_n(t-\tau) of the state that should have a delay, something like:

\frac{dy_1}{dt} = f_1(\mathbf{y}) \\
\frac{dy_2}{dt} = \rho (y_1 - y_2) \\
\frac{dy_3}{dt} = \rho (y_2 - y_3) \\
... \\
\frac{dy_n}{dt} = \rho y_{n-1} + f(\mathbf{y}) \\

that would delay the â€śavailabilityâ€ť of y_n by replacing an exponential waiting time with a gamma-distributed waiting time resulting form the sum of sequential exponential times, with mean waiting time dependent only on rate \rho and the number of compartments n (although it will increase the number of components of the state \mathbf{y}, in case you need to estimate initial conditions). This would effectively introduce a delay, although Iâ€™m not sure if and how it is related to DDEs.

If that makes sense for your problem it is formulated as an ODE without any time dependence, so it should be straightforward to do it with the Stan solver.