Stochastic differential equation and latent variables (Affine model, Vasicek model)

Modeling
,
1

If I want to fit a model that :
y=x_{1t}+x_{2t}+p_t+s_t
x_{1t},x_{2t}, p_{t}, s_t are latent variables. have distribution:
x_{it+1} \sim N([1-exp(-\kappa_i\Delta t)]\mu_i+exp(-\kappa_i)x_{it},\frac{\sigma_i^2}{2\kappa_i}[1-exp(-2\kappa_i\Delta t)]), \qquad (i=1,2)
p_t \sim N([1-exp(-\kappa_p\Delta t)]\mu_p+exp(-\kappa_p)p_{t},\frac{\sigma_p^2}{2\kappa_p}[1-exp(-2\kappa_p\Delta t)])
s_t is a stochastic volatility
The code is :

data {
  int<lower=0> T; //the number of sample
  real y1[T];     //data
}

parameters {
  //parameters for mu
  real mu1;                    //the mean for x1 
  real mu2;                    //the mean for x2
  real mup;                    //the mean for p
  real muh;                    //the mean for h
  //parameters for sigma
  real<lower=0> sigma1;        //the standard deviation(sd) for x1
  real<lower=0> sigma2;        //the sd for x2
  real<lower=0> sigmap;        //the sd for p
  real<lower=0> sigmah;        //the sd for h
  //parameters for kappa
  real<lower=0> kappa1;        
  real<lower=0> kappa2;
  real<lower=0> kappap;

  vector<lower=0>[K] sigy1; 
  real<lower=-1,upper=1> phi_h; 
  
  //parameter for intercept
  real cr;
   // latent variables
  vector[T] x1;  //latent variabe vector x1
  vector[T] x2;  //latent variabe vector x2
  vector[T] p;   //latent variabe vector p
  vector[T] h;   //latent variabe vector h
  vector[T] ep;  //latent variabe vector ep

}

transformed parameters{
  real<lower=0> sigr1;
  real<lower=0> sigr2;
  real<lower=0> sigp;
  vector[T] my1;
  sigr1=sigma1*sqrt((1-exp(-2*kappa1/52))/(2*kappa1));
  sigr2=sigma2*sqrt((1-exp(-2*kappa2/52))/(2*kappa2));
  sigp=sigmap*sqrt((1-exp(-2*kappap/52))/(2*kappap));
  
   for (t in 1:T){
    my1[t]= cr+x1[t]+x2[t]+p[t]+ep[t];
}

model{
  //set priors
  mu1 ~ normal(0,5);
  mu2 ~ normal(0,5);
  muh ~ normal(0,5);
  mup ~ normal(0,5);
  sigma1 ~ normal(0,1);
  sigma2 ~ normal(0,1);
  sigmap ~ normal(0,1);
  sigmah ~ normal(0,1);
  kappa1 ~  normal(0,10);
  kappa2 ~ normal(0,10);
  kappap ~ normal(0,10);
  phi_h ~ uniform(-1, 1);
  cr ~ normal(0,10);
  sigy1 ~ normal(0,0.1);
  //first value of latent variables
  x1[1] ~ normal(0,5);
  x2[1] ~ normal(0,5);
  p[1] ~ normal(0,5);
  h[1] ~ normal(muh, sigmah/sqrt(1-phi_h^2));
  ep[1] ~ normal(0, exp(h[1] / 2));
  for (t in 2:T){
    x1[t]~normal((1-exp(-kappa1/52))*mu1+exp(-kappa1/52)*x1[t-1],sigr1);
    x2[t]~normal((1-exp(-kappa2/52))*mu2+exp(-kappa2/52)*x2[t-1],sigr2);
    p[t] ~ normal((1-exp(-kappap/52))*mup+exp(-kappap/52)*p[t-1],sigp);
    h[t] ~ normal(muh + phi_h * (h[t - 1] -  muh), sigmah);
    ep[t] ~ normal(0, exp(h[t] / 2));
  } 
  for(t in 1:T){ 
    y1[t] ~ normal(my1[t],sigdr);
  }  
}

My questions are:
This model doesn’t have a stable solution. how can I get a stable solution, whether I need to set an initial value for SDE? Can this model be identified?

2

I don’t really understand SDEs, so I might be misrepresenting your question, but since nobody else answered, I will give it a try.

If there are multiple very different solutions, but that give you very similar likelihood, the resulting posterior will be multimodal / or have some “ridge” and this is a problem for Stan. In particular, it appears that the role of x1, x2 and p is completely symmetric and thus any permutation of all the relevant parameters will provide exactly the same result. So at the very least, you need to somehow fix the ordering, (e.g. by enforcing \kappa_1 > \kappa_2 > \kappa_p).

Additionally, I just don’t see how a single observation at each time point could let you learn something useful about 4 quantities. It appears to me, that you could introduce a new variable z_t = x_{1t} + x_{2t} + p_t and describe almost the same observed dynamics of y just in terms of z_t. It appears that the three process allow you to have somewhat more complex behaviour with regards to the speed of reversion to the overall mean, but I would be surprised

Additionally, I think ep could be completely marginalized out of the model as it is just an additional Gaussian noise that can be combined with the observation noise, e.g. by having my1[t]= cr+x1[t]+x2[t]+p[t]; and then adding the variances of the noise: y1[t] ~ normal(my1[t],sqrt(square(sigdr) + square(exp(h[t] / 2))). (this would also mean that h[1] and sigdr are not both identifiable, so you would probably want to keep just one of them in the model)

Hope that helps at least a little!