Modelling a jump component

Modeling
1

Dear All! I would like to ask a question. How to implement the following model in Stan?

I have a time-series:

y_t = \mu_t + \xi_t J_t + \epsilon_t, \qquad \epsilon_t \sim \mathcal{N}(0,\exp(h_t/2)).

I have no problem with modelling the SV part and \mu_t, but I do not understand how to incorporate a jump component. Ideally, I would like to have the following specification:

\xi_t \sim \mathcal{N}(0, \sigma_{\xi}^2) - measures the intensity of the jump.
J_t = \Phi(\lambda_0 + \lambda_x X) - determines the probability of the jump, where X is observable data that gives me information about the probability of the jump. If I use

 bernoulli(Phi_approx(lam_1  + lam_x * x));

I receive an error:

Ill-typed arguments to ‘~’ statement. No distribution ‘bernoulli’ was found with the correct signature.

Because Stan does not support discrete parameter. Could you please suggest how to translate the model in Stan model or propose a similar model that could be implemented in Stan? Below, I provide a code that I wrote so far, but as I wrote before it does not match the structure exactly

data {
  int<lower=0> T_len;           // Number of observations
  vector[T_len] y;              // Time series data
  vector[T_len] x;              // Time series data
}

parameters {
  vector[T_len] xi;                             // jump magnitude at time t
  real<lower=0> sig_y;                          // Standard deviation of the y process
  real<lower=0> sig_xi;                         // Standard deviation of the jump magnitude
  real gamma_y;                                 // Intercept of the y process
  real lam_1;                                   // Parameters for jump
  real lam_x;                                   // Parameters for jump
}

transformed parameters {
  vector[T_len] jump = Phi_approx(lam_1  + lam_x * x);          // Jump index at time t  
}

model {
  // Prior
  sig_y    ~ cauchy(0, 5);               // Prior on the standard deviation of the y process
  sig_xi   ~ cauchy(0, 5);               // Prior on the standard deviation of the jump magnitude
  gamma_y  ~ normal(0, 10);              // Prior on the intercept of the y process
  lam_1    ~ normal(0, 10);              // Prior on the intercept of the y process
  lam_x    ~ normal(0, 10);              // Prior on the intercept of the y process

  xi       ~ normal(0, sig_xi * jump);                       // Jump magnitude at time t

  // Likelihood (speed up)
  y ~ normal(gamma_y + xi, sig_y);            // Likelihood of the data given
}

Thank you very much!

Best,
Borys

2

Discrete parameters must be integrated out of Stan models.

For simplicity, let’s start with a likelihood that uses bernoulli() (incorrectly). I gather you essentially want

for (t in 1:T_len) {
  xi[t] ~ normal(0, sig_xi);
  if (bernoulli(jump[t])) {
    y[t] ~ normal(mu[t] + xi[t], sig_y);
  } else {
    y[t] ~ normal(mu[t], sig_y);
  }
}

(your code had gamma_y but I substituted mu[t] as that’s what’s shown in the equations.)

First, integrate out the jump size xi. (It’s continuous so integrating it out is not strictly required but (I expect) will improve sampling efficiency significantly). So, the likelihood is now

for (t in 1:T_len) {
  if (bernoulli(jump[t])) {
    y[t] ~ normal(mu[t], sqrt(sig_xi^2 + sig_y^2));
  } else {
    y[t] ~ normal(mu[t], sig_y);
  }
}

and the model has no parameter xi.

Integrating out discrete choices requires us to use target += syntax instead of the sampling statement. Here’s how it looks. The following code is equivalent to the previous block, just different syntax that explicitly spells out log probability density calculations.

for (t in 1:T_len) {
  if (bernoulli(jump[t])) {
    target += normal_lpdf(y[t] | mu[t], sqrt(sig_xi^2 + sig_y^2));
  } else {
    target += normal_lpdf(y[t] | mu[t], sig_y);
  }
}

Now we can handle the discrete Bernoulli choice.
Stan has a function specifically for that, log_mix.

for (t in 1:T_len) {
  target += log_mix(jump[t],
                normal_lpdf(y[t] | mu[t], sqrt(sig_xi^2 + sig_y^2)),
                normal_lpdf(y[t] | mu[t], sig_y));
}
2 Likes
3

Dear nhuurre, thank you very much for your response. It’s exactly what I looked for. I’m sorry for the late reply, since I was extremely busy and did not work on this problem