(Simplified) model
I have a system that evolves in time (discrete time), and a signal x \in \mathbb{R} can be measured at every time. The state of the system s \in \mathbb{R} has the constrain s> 0, and the signal x depends on the state as
x = as + b,
where a \in \mathbb{R} and b \sim \text{normal}(\bar{b}, \sigma_0).
The state s is only known for a few times because it is expensive to measure it directly. The model for the evolution of the state is
\Delta s_n := s_n  s_{n1} = \Delta s_{n1} + c,
where c \sim \text{normal}(\bar{c}, \sigma_1), \bar{c} < 0, and \Delta s > 0.
My task is to infer the current state s_N from x_N.
Note: the model is more complex, but I have simplified it to highlight my main questions.
Code
data {
// N is the total number of observations.
int<lower=1> N;
vector[N] x;
// Ny is the number of direct observations of the state.
// y denotes the observed states.
int<lower=1, upper=N> Ny;
array[Ny] int<lower=1, upper=N> id_y;
vector<lower=0>[Ny] y;
}
parameters {
// Parameters describing the evolution of states.
real<lower=0> s_mu, s_sigma;
// Unknown states are modelled as parameters.
vector<lower=0>[NNy] s_unknown;
// Parameters relating x and s.
real x_alpha, x_beta;
real<lower=0> x_sigma;
}
transformed parameters {
// Assign s to the known (y) and unknown (s_unknown) states.
vector<lower=0>[N] s;
{
int count_y = 1, count_s_u = 1;
for (n in 1:N) {
if (n == id_y[count_y]) {
s[n] = y[count_y];
if (count_y < Ny) {
count_y = count_y + 1;
}
}
else {
s[n] = s_unknown[count_s_u];
count_s_u = count_s_u + 1;
}
}
}
// If I remove <lower=0>, then I do not have any more problem with the sampling.
vector<lower=0>[N1] d_s;
d_s = s[2:]  s[1:N1];
}
model {
d_s[2:N1] ~ normal(s_mu + d_s[:N2], s_sigma);
x ~ normal(x_alpha * s + x_beta, x_sigma);
}
Questions
 My time series x have length of around 500, but I only know around 20 states. Is modelling the unknown s as parameters sensible, or should I try something else? Do you have some piece of advice or trick to handle this situation?
I read the page on HMM, and it is insightful, but to follow the semisupervised section I would need something like log_integral_exp, which does not seem quite right. In fact, one option is to use a HMM and to discretize the problem, but I feel my setting is more natural.

When I sample, I get errors because for some times \Delta s < 0, which is forbidden by the model.

I enclose the variables
count_
in braces because otherwise cmdstanpy complains; I am not supposed to use a parameter of typeint
. Is there other option for using local variables in the blocks.