library(nimble, warn.conflicts = FALSE)
## nimble version 0.12.2 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit https://R-nimble.org.
stochvol_installed <- require(stochvol, quietly = TRUE)

Stochastic volatility models are often used for time series of log returns of financial assets. The main idea is to model volatility (standard deviation of log returns) as an unobserved autoregressive process.

This example shows coding of a stochastic volatility model in nimble. See the particle filter and particle MCMC (PMCMC) examples for demonstrations of these methods with this model. We follow the example model and data described in the vignette (“Dealing with Stochastic Volatility in Time Series Using the R Package stochvol”) of the stochvol package (Kastner 2016, Kaster and Hosszejni 2019). We reparameterize relative to the stochvol example as noted below.

Code for this model is:

stochVolCode <- nimbleCode({
  x[1] ~ dnorm(0, sd = sigma / sqrt(1-phi*phi))
  y[1] ~ dnorm(0, sd = beta * exp(0.5 * x[1]))
  for (t in 2:T){
    x[t] ~ dnorm(phi * x[t-1], sd = sigma)
    y[t] ~ dnorm(0, sd = beta * exp(0.5 * x[t]))
  }
  phi <- 2 * phiStar - 1
  phiStar ~ dbeta(20, 1.1)
  logsigma2 ~ dgammalog(shape = 0.5, rate = 1/(2*0.1)) ## This is Omega
  sigma <- exp(0.5*logsigma2)
  mu ~ dnorm(-10, sd = 1) ## It matters whether data are converted to % or not.
  beta <- exp(0.5*mu)
})

In this model:

The probability density needed for logsigma2 is given by the following custom distributions:

dgammalog <- nimbleFunction(
  run = function(x = double(), shape = double(),
                 rate = double(),log = integer(0, default = 0)) {
    logProb <- shape * log(rate) + shape * x - rate * exp(x) - lgamma(shape)
    if(log) return(logProb)
    else return(exp(logProb))
    returnType(double())
  }
)

rgammalog <- nimbleFunction(
  run = function(n = integer(),
                 shape = double(), rate = double()) {
    xg <- rgamma(1, shape = shape, rate = rate)
    return(log(xg))
    returnType(double())
  }
)

The “r” function can be skipped if it will be not be needed, but we include it here for completeness.

Again following the stochvol vignette, we use as data exchange rates for the Euro (EUR) quoted in U.S. Dollars (USD). Here we choose to start after January 1st, 2010, and continue until the end of the time-series, 582 days after that.

message('To rebuild this example further, the stochvol package must be installed.')
## To rebuild this example further, the stochvol package must be installed.

Set up the data

data('exrates')
y <- logret(exrates$USD[exrates$date > '2010-01-01'], demean = TRUE)

Build the model

stochVolModel <- nimbleModel(code = stochVolCode,
                             constants = list(T = length(y)), data = list(y = y),
                             inits = list(mu = -10, phiStar = .99,
                                          logsigma2 = log(.004)))

Compile the model

CstochVolModel <- compileNimble(stochVolModel)

See further examples for using this model.