svsample: Markov Chain Monte Carlo (MCMC) Sampling for the Stochastic... In stochvol: Efficient Bayesian Inference for Stochastic Volatility (SV) Models

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. For similar functionality with a formula interface, see `svlm`.

Usage

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135``` ```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` may contain elements named `mu`, `phi`, `sigma`, `nu`, `rho`, `beta`, and `latent0`. 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) and Hosszejni and Kastner (2019).

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"`.

`svsample2` is deprecated.

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.

Hosszejni, D. and Kastner, G. (2019). Approaches Toward the Bayesian Estimation of the Stochastic Volatility Model with Leverage. Springer Proceedings in Mathematics & Statistics, 296, 75–83, doi: 10.1007/978-3-030-30611-3_8.

`svlm`, `svsim`, `specify_priors`
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152``` ```############### # 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(predy(preddraws), 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(predy(preddraws), 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) ```