Model Jumps in Stochastic Volatility Model

[JM]

Hi there!

I am quite new to Nimble and want to implement a stochastic volatility model with jumps.

The variable I cannot manage to implement well is defined as:

mu_t = mu_t-1 + delta*sigma*z_t

where delta is Bern(p), sigma has an inverse Gamma prior, and z_t is standard normal.

mu_t should have a normal prior.

However, I don't know how to incorporate the Bernoulli variable into mu's distribution.

Thank you so much in advance.

Best regards

JM

[Perry de Valpine]

Hi JM,

Thanks for the question. That is a fairly bare-bones question so it's a bit hard to get started with. I can point to this paper (10.18637/jss.v100.i03), which on page 20 uses a stochastic volatility model as an example. Maybe that and the other examples in the paper will help you get oriented. If you are still stuck, it would be helpful if you can try coding at least a simpler model in your problem area so that you can ask about something you want to add or modify in a more specific way. The paper is about sequential Monte Carlo methods (particle filtering), but the same model setup could be used for MCMC or other methods.

HTH

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/daab1d67-fe03-4b42-978c-27331598cae8n%40googlegroups.com.

[JM]

Hi Perry,

thanks for the reply. Below you can see the model that I want to implement.

The code that I have so far is this:

code <- nimbleCode({
  # Priors for each asset
  for (i in 1:N) {
    phi[i] ~ T(dnorm(0, sd=0.3), -0.99, 0.99) # Persistence factor
    sigma_h[i] ~ dinvgamma(1, scale=0.5)  # Std of transitory volatility
  }
  sv[1] <- 1 # Normalize first factor to 1
  for (i in 2:N) {
    sv[i] ~ dnorm(1, sd=0.5) # Scale factor for permanent component
  }

  sm[1] <- 1 # Normalize first factor to 1
  for (i in 2:N) {
    sm[i] ~ dnorm(1, sd=0.5)  # Scale factor for jumps in mean
  }
 
  # Priors for constant variables
  sigma_eta ~ dinvgamma(2, scale=2)  # Std of permanent component
  sigma_nu ~ dinvgamma(2, scale=2)  # Std of jumps in mean
  p_v ~ dbeta(2,5)  # Probability of volatility jumps
  p_m ~ dbeta(2,3)  # Probability of jumps in mean
 
  # Priors for each time step
  mu[1] ~ dnorm(-10, sd=0.2)  # Initial state of permanent volatility
  for (t in 2:T) {
    mu[t] ~ dnorm(mu[t-1], sd=delta_v[t]*sigma_eta)
  }

  for (t in 1:T) {
     u_v[t] ~ dunif(0, 1)

     u_m[t] ~ dunif(0, 1)
     z_m[t] ~ dnorm(0, 1)  # Standard normal innovations for jumps in mean

     delta_v[t] <- step(u_v[t]-(1-p_v))

     delta_m[t] <- step(u_m[t]-(1-p_m))

  }
 
  # Priors for each asset and time step
  for (i in 1:N) {
    h[1, i] ~ dnorm(0, sd=sigma_h[i]/sqrt(1-phi[i]*phi[i]))  # Initial state of transitory volatility
    for (t in 2:T) {  
      h[t, i] ~ dnorm(phi[i] * h[t-1,i], sd=sigma_h[i])  # Transitory component
    }
  }
  for (i in 1:N) {
    for (t in 1:T) {
      gamma[t, i] <- sm[i] * delta_m[t] * sigma_nu * z_m[t]  # Jump process in mean 
    }
  }
 
  # Loop over assets
  for (i in 1:N) {
    for (t in 1:T) {
      y[t, i] ~ dnorm(gamma[t, i]*exp((h[t, i] + sv[i] * mu[t]) / 2), sd = exp((h[t, i] + sv[i] * mu[t]) / 2))  # Observation equation
    }
  }
})

The parts I am not sure about are highlighted in red. Defining mu the way I did it here does not work properly, as it prevents sigma_eta from updating over iterations. So I was wondering how to properly incorporate the Bernoulli variable delta_v into the mu equation.

Further, according to the model the Bernoulli processes are sampled using auxiliary uniform variables and take on value one if the uniform variables exceed a threshold that is aligned with the parameters p_m and p_v. Therefore, I decided to use 1-p_v and 1-p_m as thresholds but I am open for other suggestions.

Lastly the observation equation contains the conditional volatility in the mean so far, which I am not sure about. Should I only use gamma as conditional mean?

Thank you very much for any suggestions!

Best regards

JM

[Perry de Valpine]

Hi JM,

I can't follow all the details (I'm not sure your code directly reflects the equations you provided, but I see it comes close) but I can try to reply.

If you want the deltas to follow a Bernoulli, you can do it directly as:

delta_v[t] ~ dbern(p_v)

Alternatively:

delta_v[t] ~ dbinom(size = 1, prob = p_v)

The declarations like

 sv[i] ~ dnorm(1, sd=0.5)

were curious to me because they are not in your equations and it seems unusual to give a prior mean of 1, but anyway that's not what you're asking about.

I guess the problem with 

mu[t] ~ dnorm(mu[t-1], sd=delta_v[t]*sigma_eta)

might occur because if delta_v[t] is 0, the log probability must be either Inf (if mu[t] == mu[t-1]) or -Inf (otherwise). 

MCMC updates of the mu[t] values will never make them perfectly equal, so you'll never be able to have delta_v[t] == 0 accepted in MCMC. (Also I don't see eta_t represented in your code.)

You might be able to set this up by adding the second term of equation (3) to mu[t-1] and using the reversible jump (see configureRJ) with delta_v[t] as the indicator variable. I'm not completely sure. That does create a potentially inefficient way to write a state-space model in nimble code because then the sequence of states would be deterministic sums of stochastic nodes, such that updating at time t requires would calculation of all subsequent times. Another more involved option would be to write a joint MCMC sampler for delta_v[t] and mu[t], if you are comfortable creating MCMC samplers. The User Manual covers the software side of this in section 15.5.

The questions about how to use an auxiliary variable (if you don't want to use dbern or dbinom directly) or the correct observation model seem to be more about the statistics and not the software, and I don't know these models well enough to respond. Someone else might chime in, but generally our support on this list is on the software mechanics. Let me know if I'm not seeing the questions clearly.

Perry

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/30159cf6-15a3-4c7c-bd88-f5a41b30f02en%40googlegroups.com.