stochvol: svsample – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

svsample

Markov Chain Monte Carlo (MCMC) Sampling for the Stochastic Volatility (SV) Model

Description

svsample simulates from the joint posterior distribution of the SV parameters mu, phi, sigma (and potentially nu and rho), along with the latent log-volatilities h_0,...,h_n and returns the MCMC draws. If a design matrix is provided, simple Bayesian regression can also be conducted.

Usage

svsample(
  y,
  draws = 10000,
  burnin = 1000,
  designmatrix = NA,
  priormu = c(0, 100),
  priorphi = c(5, 1.5),
  priorsigma = 1,
  priornu = 0,
  priorrho = NA,
  priorbeta = c(0, 10000),
  priorlatent0 = "stationary",
  priorspec = NULL,
  thin = 1,
  thinpara = thin,
  thinlatent = thin,
  keeptime = "all",
  quiet = FALSE,
  startpara = NULL,
  startlatent = NULL,
  parallel = c("no", "multicore", "snow"),
  n_cpus = 1L,
  cl = NULL,
  n_chains = 1L,
  print_progress = "automatic",
  expert = NULL,
  ...
)

svtsample(
  y,
  draws = 10000,
  burnin = 1000,
  designmatrix = NA,
  priormu = c(0, 100),
  priorphi = c(5, 1.5),
  priorsigma = 1,
  priornu = 0.1,
  priorrho = NA,
  priorbeta = c(0, 10000),
  priorlatent0 = "stationary",
  priorspec = NULL,
  thin = 1,
  thinpara = thin,
  thinlatent = thin,
  keeptime = "all",
  quiet = FALSE,
  startpara = NULL,
  startlatent = NULL,
  parallel = c("no", "multicore", "snow"),
  n_cpus = 1L,
  cl = NULL,
  n_chains = 1L,
  print_progress = "automatic",
  expert = NULL,
  ...
)

svlsample(
  y,
  draws = 20000,
  burnin = 2000,
  designmatrix = NA,
  priormu = c(0, 100),
  priorphi = c(5, 1.5),
  priorsigma = 1,
  priornu = 0,
  priorrho = c(4, 4),
  priorbeta = c(0, 10000),
  priorlatent0 = "stationary",
  priorspec = NULL,
  thin = 1,
  thinpara = thin,
  thinlatent = thin,
  keeptime = "all",
  quiet = FALSE,
  startpara = NULL,
  startlatent = NULL,
  parallel = c("no", "multicore", "snow"),
  n_cpus = 1L,
  cl = NULL,
  n_chains = 1L,
  print_progress = "automatic",
  expert = NULL,
  ...
)

svtlsample(
  y,
  draws = 20000,
  burnin = 2000,
  designmatrix = NA,
  priormu = c(0, 100),
  priorphi = c(5, 1.5),
  priorsigma = 1,
  priornu = 0.1,
  priorrho = c(4, 4),
  priorbeta = c(0, 10000),
  priorlatent0 = "stationary",
  priorspec = NULL,
  thin = 1,
  thinpara = thin,
  thinlatent = thin,
  keeptime = "all",
  quiet = FALSE,
  startpara = NULL,
  startlatent = NULL,
  parallel = c("no", "multicore", "snow"),
  n_cpus = 1L,
  cl = NULL,
  n_chains = 1L,
  print_progress = "automatic",
  expert = NULL,
  ...
)

svsample2(
  y,
  draws = 10000,
  burnin = 1000,
  designmatrix = NA,
  priormu = c(0, 100),
  priorphi = c(5, 1.5),
  priorsigma = 1,
  priornu = 0,
  priorrho = NA,
  priorbeta = c(0, 10000),
  priorlatent0 = "stationary",
  thinpara = 1,
  thinlatent = 1,
  keeptime = "all",
  quiet = FALSE,
  startpara = NULL,
  startlatent = NULL
)

Arguments

`y`	numeric vector containing the data (usually log-returns), which must not contain zeros. Alternatively, `y` can be an `svsim` object. In this case, the returns will be extracted and a message is signalled.
`draws`	single number greater or equal to 1, indicating the number of draws after burn-in (see below). Will be automatically coerced to integer. The default value is 10000.
`burnin`	single number greater or equal to 0, indicating the number of draws discarded as burn-in. Will be automatically coerced to integer. The default value is 1000.
`designmatrix`	regression design matrix for modeling the mean. Must have `length(y)` rows. Alternatively, `designmatrix` may be a string of the form `"arX"`, where `X` is a nonnegative integer. To fit a constant mean model, use `designmatrix = "ar0"` (which is equivalent to `designmatrix = matrix(1, nrow = length(y))`). To fit an AR(1) model, use `designmatrix = "ar1"`, and so on. If some elements of `designmatrix` are `NA`, the mean is fixed to zero (pre-1.2.0 behavior of stochvol).
`priormu`	numeric vector of length 2, indicating mean and standard deviation for the Gaussian prior distribution of the parameter `mu`, the level of the log-volatility. The default value is `c(0, 100)`, which constitutes a practically uninformative prior for common exchange rate datasets, stock returns and the like.
`priorphi`	numeric vector of length 2, indicating the shape parameters for the Beta prior distribution of the transformed parameter `(phi + 1) / 2`, where `phi` denotes the persistence of the log-volatility. The default value is `c(5, 1.5)`, which constitutes a prior that puts some belief in a persistent log-volatility but also encompasses the region where `phi` is around 0.
`priorsigma`	single positive real number, which stands for the scaling of the transformed parameter `sigma^2`, where `sigma` denotes the volatility of log-volatility. More precisely, `sigma^2 ~ priorsigma * chisq(df = 1)`. The default value is `1`, which constitutes a reasonably vague prior for many common exchange rate datasets, stock returns and the like.
`priornu`	single non-negative number, indicating the rate parameter for the exponential prior distribution of the parameter; can be `Inf` `nu`, the degrees-of-freedom parameter of the conditional innovations t-distribution. The default value is `0`, fixing the degrees-of-freedom to infinity. This corresponds to conditional standard normal innovations, the pre-1.1.0 behavior of stochvol.
`priorrho`	either `NA` for the no-leverage case or a numeric vector of length 2 that specify the beta prior distribution for `(rho+1)/2`
`priorbeta`	numeric vector of length 2, indicating the mean and standard deviation of the Gaussian prior for the regression parameters. The default value is `c(0, 10000)`, which constitutes a very vague prior for many common datasets. Not used if `designmatrix` is `NA`.
`priorlatent0`	either a single non-negative number or the string `'stationary'` (the default, also the behavior before version 1.3.0). When `priorlatent0` is equal to `'stationary'`, the stationary distribution of the latent AR(1)-process is used as the prior for the initial log-volatility `h_0`. When `priorlatent0` is equal to a number B, we have h_0 \sim N(μ, Bσ^2) a priori.
`priorspec`	in case one needs different prior distributions than the ones specified by `priormu`, `...`, `priorrho`, a `priorspec` object can be supplied here. A smart constructor for this usecase is specify_priors.
`thin`	single number greater or equal to 1, coercible to integer. Every `thinpara`th parameter and latent draw is kept and returned. The default value is 1, corresponding to no thinning of the parameter draws i.e. every draw is stored.
`thinpara`	single number greater or equal to 1, coercible to integer. Every `thinpara`th parameter draw is kept and returned. The default value is `thin`.
`thinlatent`	single number greater or equal to 1, coercible to integer. Every `thinlatent`th latent variable draw is kept and returned. The default value is `thin`
`keeptime`	Either 'all' (the default) or 'last'. Indicates which latent volatility draws should be stored.
`quiet`	logical value indicating whether the progress bar and other informative output during sampling should be omitted. The default value is `FALSE`, implying verbose output.
`startpara`	optional named list, containing the starting values for the parameter draws. If supplied, `startpara` must contain three elements named `mu`, `phi`, and `sigma`, where `mu` is an arbitrary numerical value, `phi` is a real number between `-1` and `1`, and `sigma` is a positive real number. Moreover, if `priornu` is not `0`, `startpara` must also contain an element named `nu` (the degrees of freedom parameter for the t-innovations). The default value is equal to the prior mean. In case of parallel execution with `cl` provided, `startpara` can be a list of named lists that initialize the parallel chains.
`startlatent`	optional vector of length `length(y)`, containing the starting values for the latent log-volatility draws. The default value is `rep(-10, length(y))`. In case of parallel execution with `cl` provided, `startlatent` can be a list of named lists that initialize the parallel chains.
`parallel`	optional one of `"no"` (default), `"multicore"`, or `"snow"`, indicating what type of parallellism is to be applied. Option `"multicore"` is not available on Windows.
`n_cpus`	optional positive integer, the number of CPUs to be used in case of parallel computations. Defaults to `1L`. Ignored if parameter `cl` is supplied and `parallel != "snow"`.
`cl`	optional so-called SNOW cluster object as implemented in package `parallel`. Ignored unless `parallel == "snow"`.
`n_chains`	optional positive integer specifying the number of independent MCMC chains
`print_progress`	optional one of `"automatic"`, `"progressbar"`, or `"iteration"`, controls the output. Ignored if `quiet` is `TRUE`.
`expert`	optional named list of expert parameters. For most applications, the default values probably work best. Interested users are referred to the literature provided in the References section. If `expert` is provided, it may contain the following named elements: interweaveLogical value. If `TRUE` (the default), then ancillarity-sufficiency interweaving strategy (ASIS) is applied to improve on the sampling efficiency for the parameters. Otherwise one parameterization is used. correct_model_misspecificationLogical value. If `FALSE` (the default), then auxiliary mixture sampling is used to sample the latent states. If `TRUE`, extra computations are made to correct for model misspecification either ex-post by reweighting or on-line using a Metropolis-Hastings step.
`...`	Any extra arguments will be forwarded to `updatesummary`, controlling the type of statistics calculated for the posterior draws.

Details

Functions svtsample, svlsample, and svtlsample are wrappers around svsample with convenient default values for the SV model with t-errors, leverage, and both t-errors and leverage, respectively.

For details concerning the algorithm please see the paper by Kastner and Frühwirth-Schnatter (2014).

Value

The value returned is a list object of class svdraws holding

`para`	`mcmc.list` object containing the parameter draws from the posterior distribution.
`latent`	`mcmc.list` object containing the latent instantaneous log-volatility draws from the posterior distribution.
`latent0`	`mcmc.list` object containing the latent initial log-volatility draws from the posterior distribution.
`tau`	`mcmc.list` object containing the latent variance inflation factors for the sampler with conditional t-innovations (optional).
`beta`	`mcmc.list` object containing the regression coefficient draws from the posterior distribution (optional).
`y`	the left hand side of the observation equation, usually the argument `y`. In case of an AR(`k`) specification, the first `k` elements are removed.
`runtime`	`proc_time` object containing the run time of the sampler.
`priors`	a `priorspec` object containing the parameter values of the prior distributions for `mu`, `phi`, `sigma`, `nu`, `rho`, and `beta`s, and the variance of specification for `latent0`.
`thinning`	`list` containing the thinning parameters, i.e. the arguments `thinpara`, `thinlatent` and `keeptime`.
`summary`	`list` containing a collection of summary statistics of the posterior draws for `para`, `latent`, and `latent0`.
`meanmodel`	`character` containing information about how `designmatrix` was employed.

To display the output, use print, summary and plot. The print method simply prints the posterior draws (which is very likely a lot of output); the summary method displays the summary statistics currently stored in the object; the plot method plot.svdraws gives a graphical overview of the posterior distribution by calling volplot, traceplot and densplot and displaying the results on a single page.

Note

If y contains zeros, you might want to consider de-meaning your returns or use designmatrix = "ar0".

References

Kastner, G. and Frühwirth-Schnatter, S. (2014). Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models. Computational Statistics & Data Analysis, 76, 408–423, doi: 10.1016/j.csda.2013.01.002.

Examples

###############
# Full examples
###############

# Example 1
## Simulate a short and highly persistent SV process 
sim <- svsim(100, mu = -10, phi = 0.99, sigma = 0.2)

## Obtain 5000 draws from the sampler (that's not a lot)
draws <-
  svsample(sim, draws = 5000, burnin = 100,
           priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2)

## Check out the results
summary(draws)
plot(draws)



# Example 2
## Simulate an asymmetric and conditionally heavy-tailed SV process
sim <- svsim(150, mu = -10, phi = 0.96, sigma = 0.3, nu = 10, rho = -0.3)

## Obtain 10000 draws from the sampler
## Use more advanced prior settings
## Run two parallel MCMC chains
advanced_draws <-
  svsample(sim, draws = 10000, burnin = 5000,
           priorspec = specify_priors(mu = sv_normal(-10, 1),
                                      sigma2 = sv_gamma(0.5, 2),
                                      rho = sv_beta(4, 4),
                                      nu = sv_constant(5)),
           parallel = "snow", n_chains = 2, n_cpus = 2)

## Check out the results
summary(advanced_draws)
plot(advanced_draws)


# Example 3
## AR(1) structure for the mean
data(exrates)
len <- 3000
ahead <- 100
y <- head(exrates$USD, len)

## Fit AR(1)-SVL model to EUR-USD exchange rates
res <- svsample(y, designmatrix = "ar1")

## Use predict.svdraws to obtain predictive distributions
preddraws <- predict(res, steps = ahead)

## Calculate predictive quantiles
predquants <- apply(preddraws$y[[1]], 2, quantile, c(.1, .5, .9))

## Visualize
expost <- tail(head(exrates$USD, len+ahead), ahead)
ts.plot(y, xlim = c(length(y)-4*ahead, length(y)+ahead),
	       ylim = range(c(predquants, expost, tail(y, 4*ahead))))
for (i in 1:3) {
  lines((length(y)+1):(length(y)+ahead), predquants[i,],
        col = 3, lty = c(2, 1, 2)[i])
}
lines((length(y)+1):(length(y)+ahead), expost,
      col = 2)


# Example 4
## Predicting USD based on JPY and GBP in the mean
data(exrates)
len <- 3000
ahead <- 30
## Calculate log-returns
logreturns <- apply(exrates[, c("USD", "JPY", "GBP")], 2,
                    function (x) diff(log(x)))
logretUSD <- logreturns[2:(len+1), "USD"]
regressors <- cbind(1, as.matrix(logreturns[1:len, ]))  # lagged by 1 day

## Fit SV model to EUR-USD exchange rates
res <- svsample(logretUSD, designmatrix = regressors)

## Use predict.svdraws to obtain predictive distributions
predregressors <- cbind(1, as.matrix(logreturns[(len+1):(len+ahead), ]))
preddraws <- predict(res, steps = ahead,
                     newdata = predregressors)
predprice <- exrates[len+2, "USD"] * exp(t(apply(preddraws$y[[1]], 1, cumsum)))

## Calculate predictive quantiles
predquants <- apply(predprice, 2, quantile, c(.1, .5, .9))

## Visualize
priceUSD <- exrates[3:(len+2), "USD"]
expost <- exrates[(len+3):(len+ahead+2), "USD"]
ts.plot(priceUSD, xlim = c(len-4*ahead, len+ahead+1),
	       ylim = range(c(expost, predquants, tail(priceUSD, 4*ahead))))
for (i in 1:3) {
  lines(len:(len+ahead), c(tail(priceUSD, 1), predquants[i,]),
        col = 3, lty = c(2, 1, 2)[i])
}
lines(len:(len+ahead), c(tail(priceUSD, 1), expost),
      col = 2)


########################
# Further short examples
########################

y <- svsim(50, nu = 10, rho = -0.1)$y

# Supply initial values
res <- svsample(y,
                startpara = list(mu = -10, sigma = 1))


# Supply initial values for parallel chains
res <- svsample(y,
                startpara = list(list(mu = -10, sigma = 1),
                                 list(mu = -11, sigma = .1, phi = 0.9),
                                 list(mu = -9, sigma = .3, phi = 0.7)),
                parallel = "snow", n_chains = 3, n_cpus = 2)

# Parallel chains with with a pre-defined cluster object
cl <- parallel::makeCluster(2, "PSOCK", outfile = NULL)  # print to console
res <- svsample(y,
                parallel = "snow", n_chains = 3, cl = cl)
parallel::stopCluster(cl)


# Turn on correction for model misspecification
## Since the approximate model is fast and it is working very
##   well in practice, this is turned off by default
res <- svsample(y,
                expert = list(correct_model_misspecification = TRUE))
print(res)

## Not run: 
# Parallel multicore chains (not available on Windows)
res <- svsample(y, draws = 30000, burnin = 10000,
                parallel = "multicore", n_chains = 3, n_cpus = 2)

# Plot using a color palette
palette(rainbow(coda::nchain(para(res, "all"))))
plot(res)

# Use functionality from package 'coda'
## E.g. Geweke's convergence diagnostics
coda::geweke.diag(para(res, "all")[, c("mu", "phi", "sigma")])

# Use functionality from package 'bayesplot'
bayesplot::mcmc_pairs(res, pars = c("sigma", "mu", "phi", "h_0", "h_15"))

## End(Not run)

stochvol

Efficient Bayesian Inference for Stochastic Volatility (SV) Models

v3.0.4

GPL (>= 2)

Authors

Darjus Hosszejni [aut, cre] (<https://orcid.org/0000-0002-3803-691X>), Gregor Kastner [aut] (<https://orcid.org/0000-0002-8237-8271>)

Initial release