, Gregor Kastner [aut]
| Title: | Efficient Bayesian Inference for Stochastic Volatility (SV) Models |
|---|---|
| Description: | Efficient algorithms for fully Bayesian estimation of stochastic volatility (SV) models with and without asymmetry (leverage) via Markov chain Monte Carlo (MCMC) methods. Methodological details are given in Kastner and Frühwirth-Schnatter (2014) <doi:10.1016/j.csda.2013.01.002> and Hosszejni and Kastner (2019) <doi:10.1007/978-3-030-30611-3_8>; the most common use cases are described in Hosszejni and Kastner (2021) <doi:10.18637/jss.v100.i12> and Kastner (2016) <doi:10.18637/jss.v069.i05> and the package examples. |
| Authors: | Darjus Hosszejni [aut, cre]
, Gregor Kastner [aut]
|
| Maintainer: | Darjus Hosszejni <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 3.2.5 |
| Built: | 2024-11-28 06:57:00 UTC |
| Source: | CRAN |
This package provides an efficient algorithm for fully Bayesian estimation of stochastic volatility (SV) models via Markov chain Monte Carlo (MCMC) methods. Methodological details are given in Kastner and Frühwirth-Schnatter (2014); the most common use cases are described in Kastner (2016). Recently, the package has been extended to allow for the leverage effect.
Bayesian inference for stochastic volatility models using MCMC methods highly depends on actual parameter values in terms of sampling efficiency. While draws from the posterior utilizing the standard centered parameterization break down when the volatility of volatility parameter in the latent state equation is small, non-centered versions of the model show deficiencies for highly persistent latent variable series. The novel approach of ancillarity-sufficiency interweaving (Yu and Meng, 2011) has recently been shown to aid in overcoming these issues for a broad class of multilevel models. This package provides software for “combining best of different worlds” which allows for inference for parameter constellations that have previously been infeasible to estimate without the need to select a particular parameterization beforehand.
This package is currently in active development. Your comments, suggestions and requests are warmly welcome!
Gregor Kastner [email protected], Darjus Hosszejni [email protected]
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.
Kastner, G. (2016). Dealing with Stochastic Volatility in Time Series Using the R Package stochvol. Journal of Statistical Software, 69(5), 1–30, doi:10.18637/jss.v069.i05.
Yu, Y. and Meng, X.-L. (2011). To Center or Not to Center: That is Not the Question—An Ancillarity-Suffiency Interweaving Strategy (ASIS) for Boosting MCMC Efficiency. Journal of Computational and Graphical Statistics, 20(3), 531–570, doi:10.1198/jcgs.2011.203main.
Useful links:
Report bugs at https://github.com/gregorkastner/stochvol/issues
## Simulate a highly persistent SV process sim <- svsim(500, mu = -10, phi = 0.99, sigma = 0.2) ## Obtain 4000 draws from the sampler (that's too few!) draws <- svsample(sim$y, draws = 4000, burnin = 100, priormu = c(-10, 1), priorphi = c(20, 1.2), priorsigma = 0.2) ## Predict 20 days ahead fore <- predict(draws, 20) ## plot the results plot(draws, forecast = fore) ## Not run: ## Simulate an SV process with leverage sim <- svsim(500, mu = -10, phi = 0.95, sigma = 0.2, rho=-0.5) ## Obtain 8000 draws from the sampler (that's too little!) draws <- svsample(sim$y, draws = 4000, burnin = 3000, priormu = c(-10, 1), priorphi = c(20, 1.2), priorsigma = 0.2, priorrho = c(1, 1)) ## Predict 20 days ahead fore <- predict(draws, 20) ## plot the results plot(draws, forecast = fore) ## End(Not run)## Simulate a highly persistent SV process sim <- svsim(500, mu = -10, phi = 0.99, sigma = 0.2) ## Obtain 4000 draws from the sampler (that's too few!) draws <- svsample(sim$y, draws = 4000, burnin = 100, priormu = c(-10, 1), priorphi = c(20, 1.2), priorsigma = 0.2) ## Predict 20 days ahead fore <- predict(draws, 20) ## plot the results plot(draws, forecast = fore) ## Not run: ## Simulate an SV process with leverage sim <- svsim(500, mu = -10, phi = 0.95, sigma = 0.2, rho=-0.5) ## Obtain 8000 draws from the sampler (that's too little!) draws <- svsample(sim$y, draws = 4000, burnin = 3000, priormu = c(-10, 1), priorphi = c(20, 1.2), priorsigma = 0.2, priorrho = c(1, 1)) ## Predict 20 days ahead fore <- predict(draws, 20) ## plot the results plot(draws, forecast = fore) ## End(Not run)
The data set contains the daily bilateral prices of one Euro in 23 currencies from January 3, 2000, until April 4, 2012. Conversions to New Turkish Lira and Fourth Romanian Leu have been incorporated.
ECB Statistical Data Warehouse (https://sdw.ecb.europa.eu)
## Not run: data(exrates) dat <- logret(exrates$USD, demean = TRUE) ## de-meaned log-returns res <- svsample(dat) ## run MCMC sampler plot(res, forecast = 100) ## display results ## End(Not run)## Not run: data(exrates) dat <- logret(exrates$USD, demean = TRUE) ## de-meaned log-returns res <- svsample(dat) ## run MCMC sampler plot(res, forecast = 100) ## display results ## End(Not run)
Some simple extractors returning the corresponding element of an
svdraws and svpredict object.
para(x, chain = "concatenated") latent0(x, chain = "concatenated") latent(x, chain = "concatenated") vola(x, chain = "concatenated") svbeta(x, chain = "concatenated") svtau(x, chain = "concatenated") priors(x) thinning(x) runtime(x) sampled_parameters(x) predy(y, chain = "concatenated") predlatent(y, chain = "concatenated") predvola(y, chain = "concatenated")para(x, chain = "concatenated") latent0(x, chain = "concatenated") latent(x, chain = "concatenated") vola(x, chain = "concatenated") svbeta(x, chain = "concatenated") svtau(x, chain = "concatenated") priors(x) thinning(x) runtime(x) sampled_parameters(x) predy(y, chain = "concatenated") predlatent(y, chain = "concatenated") predvola(y, chain = "concatenated")
x |
|
chain |
optional either a positive integer or the string
|
y |
|
The return value depends on the actual funtion.
para(x, chain = "concatenated")
|
extracts the parameter draws. |
latent(x, chain = "concatenated")
|
extracts the latent contemporaneous log-volatility draws. |
latent0(x, chain = "concatenated")
|
extracts the latent initial log-volatility draws. |
svbeta(x, chain = "concatenated")
|
extracts the linear regression coefficient draws. |
svtau(x, chain = "concatenated")
|
extracts the tau draws. |
vola(x, chain = "concatenated")
|
extracts standard deviation draws. |
priors(x) |
extracts the prior
parameters used and returns them in a |
thinning(x) |
extracts the thinning parameters used and returns them in
a |
runtime(x) |
extracts the runtime and returns it as a
|
sampled_parameters(x) |
returns the names of time independent model
parameters that were actually sampled by |
predlatent(y, chain = "concatenated")
|
extracts the predicted latent contemporaneous log-volatility draws. |
predvola(y, chain = "concatenated")
|
extracts predicted standard deviation draws. |
predy(y, chain = "concatenated")
|
extracts the predicted observation draws. |
Functions that have input parameter chain return
an mcmc.list object if chain=="all" and
return an mcmc object otherwise. If chain is
an integer, then the specified chain is selected from
all chains. If chain is "concatenated",
then all chains are merged into one mcmc object.
specify_priors, svsample, predict.svdraws
# Simulate data sim <- svsim(150) # Draw from vanilla SV draws <- svsample(sim, draws = 2000) ## Summarize all estimated parameter draws as a merged mcmc object summary(para(draws)[, sampled_parameters(draws)]) ## Extract the draws as an mcmc.list object params <- para(draws, chain = "all")[, sampled_parameters(draws)] options(max.print = 100) ## Further short examples summary(latent0(draws)) summary(latent(draws)) summary(vola(draws)) sampled_parameters(draws) priors(draws) # Draw 3 independent chains from heavy-tailed and asymmetric SV with AR(2) structure draws <- svsample(sim, draws = 20000, burnin = 3000, designmatrix = "ar2", priornu = 0.1, priorrho = c(4, 4), n_chains = 3) ## Extract beta draws from the second chain svbeta(draws, chain = 2) ## ... tau draws from all chains merged/concatenated together svtau(draws) ## Create a new svdraws object from the first and third chain second_chain_excluded <- draws[c(1, 3)] # Draw from the predictive distribution pred <- predict(draws, steps = 2) ## Extract the predicted observations as an mcmc.list object predicted_y <- predy(pred, chain = "all") ## ... the predicted standard deviations from the second chain predicted_sd <- predvola(pred, chain = 2) ## Create a new svpredict object from the first and third chain second_chain_excluded <- pred[c(1, 3)]# Simulate data sim <- svsim(150) # Draw from vanilla SV draws <- svsample(sim, draws = 2000) ## Summarize all estimated parameter draws as a merged mcmc object summary(para(draws)[, sampled_parameters(draws)]) ## Extract the draws as an mcmc.list object params <- para(draws, chain = "all")[, sampled_parameters(draws)] options(max.print = 100) ## Further short examples summary(latent0(draws)) summary(latent(draws)) summary(vola(draws)) sampled_parameters(draws) priors(draws) # Draw 3 independent chains from heavy-tailed and asymmetric SV with AR(2) structure draws <- svsample(sim, draws = 20000, burnin = 3000, designmatrix = "ar2", priornu = 0.1, priorrho = c(4, 4), n_chains = 3) ## Extract beta draws from the second chain svbeta(draws, chain = 2) ## ... tau draws from all chains merged/concatenated together svtau(draws) ## Create a new svdraws object from the first and third chain second_chain_excluded <- draws[c(1, 3)] # Draw from the predictive distribution pred <- predict(draws, steps = 2) ## Extract the predicted observations as an mcmc.list object predicted_y <- predy(pred, chain = "all") ## ... the predicted standard deviations from the second chain predicted_sd <- predvola(pred, chain = 2) ## Create a new svpredict object from the first and third chain second_chain_excluded <- pred[c(1, 3)]
These functions define meaningful expert settings for
argument expert of svsample and its derivatives.
The result of get_default_fast_sv should be provided as
expert$fast_sv and get_default_general_sv as expert$general_sv
when relevant.
get_default_fast_sv() get_default_general_sv(priorspec) default_fast_svget_default_fast_sv() get_default_general_sv(priorspec) default_fast_sv
priorspec |
a |
An object of class list of length 11.
default_fast_sv is deprecated.
svsample, specify_priors, svsample_roll, svsample_fast_cpp, svsample_general_cpp
logret computes the log returns of a time
series, with optional de-meaning and/or standardization.
logret(dat, demean = FALSE, standardize = FALSE, ...) ## Default S3 method: logret(dat, demean = FALSE, standardize = FALSE, ...)logret(dat, demean = FALSE, standardize = FALSE, ...) ## Default S3 method: logret(dat, demean = FALSE, standardize = FALSE, ...)
dat |
The raw data. |
demean |
Logical value indicating whether the data should be de-meaned. |
standardize |
Logical value indicating whether the data should be standardized (in the sense that each component series has an empirical variance equal to one). |
... |
Ignored. |
Log returns of the (de-meaned / standardized) data.
logret(default): Log returns of vectors
Displays a plot of the density estimate for the posterior distribution of
the parameters mu, phi, sigma (and potentially
nu or rho), computed by the density function.
paradensplot( x, showobs = TRUE, showprior = TRUE, showxlab = TRUE, mar = c(1.9, 1.9, 1.9, 0.5), mgp = c(2, 0.6, 0), simobj = NULL, ... )paradensplot( x, showobs = TRUE, showprior = TRUE, showxlab = TRUE, mar = c(1.9, 1.9, 1.9, 0.5), mgp = c(2, 0.6, 0), simobj = NULL, ... )
x |
|
showobs |
logical value, indicating whether the observations should be
displayed along the x-axis. If many draws have been obtained, the default
( |
showprior |
logical value, indicating whether the prior distribution
should be displayed. The default value is |
showxlab |
logical value, indicating whether the x-axis should be
labelled with the number of iterations and the bandwith obtained from
|
mar |
numerical vector of length 4, indicating the plot margins. See
|
mgp |
numerical vector of length 3, indicating the axis and label
positions. See |
simobj |
object of class |
... |
further arguments are passed on to the invoked |
paradensplot is modeled after densplot in the
coda package, with some modifications for parameters that have
(half-)bounded support.
Called for its side effects. Returns argument x invisibly.
You can call this function directly, but it is more commonly called by
the plot.svdraws method.
Other plotting:
paratraceplot(),
paratraceplot.svdraws(),
plot.svdraws(),
plot.svpredict(),
volplot()
Generic function for plotting iterations vs. sampled parameter values.
A detailed help for the method implemented in stochvol can be found in
paratraceplot.svdraws.
paratraceplot(x, ...)paratraceplot(x, ...)
x |
An object used to select a method. |
... |
Further arguments passed to or from other methods. |
Called for its side effects. Returns argument x invisibly.
Other plotting:
paradensplot(),
paratraceplot.svdraws(),
plot.svdraws(),
plot.svpredict(),
volplot()
Displays a plot of iterations vs. sampled values the parameters mu,
phi, sigma (and potentially nu or rho), with a separate plot
per variable.
## S3 method for class 'svdraws' paratraceplot( x, mar = c(1.9, 1.9, 1.9, 0.5), mgp = c(2, 0.6, 0), simobj = NULL, ... )## S3 method for class 'svdraws' paratraceplot( x, mar = c(1.9, 1.9, 1.9, 0.5), mgp = c(2, 0.6, 0), simobj = NULL, ... )
x |
|
mar |
numerical vector of length 4, indicating the plot margins. See
|
mgp |
numerical vector of length 3, indicating the axis and label
positions. See |
simobj |
object of class |
... |
further arguments are passed on to the invoked |
paratraceplot is modeled after traceplot in the
coda package, with very minor modifications.
Called for its side effects. Returns argument x invisibly.
You can call this function directly, but it is more commonly called by
the plot.svdraws method.
Other plotting:
paradensplot(),
paratraceplot(),
plot.svdraws(),
plot.svpredict(),
volplot()
plot.svdraws and plot.svldraws generate some plots visualizing the posterior
distribution and can also be used to display predictive distributions of
future volatilities.
## S3 method for class 'svdraws' plot( x, forecast = NULL, dates = NULL, show0 = FALSE, showobs = TRUE, showprior = TRUE, forecastlty = NULL, tcl = -0.4, mar = c(1.9, 1.9, 1.7, 0.5), mgp = c(2, 0.6, 0), simobj = NULL, newdata = NULL, ... )## S3 method for class 'svdraws' plot( x, forecast = NULL, dates = NULL, show0 = FALSE, showobs = TRUE, showprior = TRUE, forecastlty = NULL, tcl = -0.4, mar = c(1.9, 1.9, 1.7, 0.5), mgp = c(2, 0.6, 0), simobj = NULL, newdata = NULL, ... )
x |
|
forecast |
nonnegative integer or object of class |
dates |
vector of length |
show0 |
logical value, indicating whether the initial volatility
|
showobs |
logical value, indicating whether the observations should be
displayed along the x-axis. If many draws have been obtained, the default
( |
showprior |
logical value, indicating whether the prior distribution
should be displayed. The default value is |
forecastlty |
vector of line type values (see
|
tcl |
The length of tick marks as a fraction of the height of a line of
text. See |
mar |
numerical vector of length 4, indicating the plot margins. See
|
mgp |
numerical vector of length 3, indicating the axis and label
positions. See |
simobj |
object of class |
newdata |
corresponds to parameter |
... |
further arguments are passed on to the invoked plotting functions. |
These functions set up the page layout and call volplot,
paratraceplot and paradensplot.
Called for its side effects. Returns argument x invisibly.
In case you want different quantiles to be plotted, use
updatesummary on the svdraws object first. An example
of doing so is given in the Examples section.
Gregor Kastner [email protected]
updatesummary, predict.svdraws
Other plotting:
paradensplot(),
paratraceplot(),
paratraceplot.svdraws(),
plot.svpredict(),
volplot()
## 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$y, draws = 5000, burnin = 1000, priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2) ## Plot the latent volatilities and some forecasts plot(draws, forecast = 10) ## Re-plot with different quantiles newquants <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99) draws <- updatesummary(draws, quantiles = newquants) plot(draws, forecast = 20, showobs = FALSE, forecastlty = 3, showprior = FALSE)## 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$y, draws = 5000, burnin = 1000, priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2) ## Plot the latent volatilities and some forecasts plot(draws, forecast = 10) ## Re-plot with different quantiles newquants <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99) draws <- updatesummary(draws, quantiles = newquants) plot(draws, forecast = 20, showobs = FALSE, forecastlty = 3, showprior = FALSE)
plot.svpredict and plot.svlpredict generate some plots
visualizing the posterior predictive distribution of future volatilites and
future observations.
## S3 method for class 'svpredict' plot(x, quantiles = c(0.05, 0.25, 0.5, 0.75, 0.95), ...)## S3 method for class 'svpredict' plot(x, quantiles = c(0.05, 0.25, 0.5, 0.75, 0.95), ...)
x |
|
quantiles |
Which quantiles to plot? Defaults to
|
... |
further arguments are passed on to the invoked
|
Called for its side effects. Returns argument x invisibly.
Note that svpredict or svlpredict objects can also be
used within plot.svdraws for a possibly more useful
visualization. See the examples in predict.svdraws and
those below for use cases.
Other plotting:
paradensplot(),
paratraceplot(),
paratraceplot.svdraws(),
plot.svdraws(),
volplot()
Other plotting:
paradensplot(),
paratraceplot(),
paratraceplot.svdraws(),
plot.svdraws(),
volplot()
## Simulate a short and highly persistent SV process sim <- svsim(100, mu = -10, phi = 0.99, sigma = 0.1) ## Obtain 5000 draws from the sampler (that's not a lot) draws <- svsample(sim$y, draws = 5000, burnin = 1000) ## Predict 10 steps ahead pred <- predict(draws, 10) ## Visualize the predicted distributions plot(pred) ## Plot the latent volatilities and some forecasts plot(draws, forecast = pred)## Simulate a short and highly persistent SV process sim <- svsim(100, mu = -10, phi = 0.99, sigma = 0.1) ## Obtain 5000 draws from the sampler (that's not a lot) draws <- svsample(sim$y, draws = 5000, burnin = 1000) ## Predict 10 steps ahead pred <- predict(draws, 10) ## Visualize the predicted distributions plot(pred) ## Plot the latent volatilities and some forecasts plot(draws, forecast = pred)
Simulates draws from the predictive density of the returns and the latent log-volatility process. The same mean model is used for prediction as was used for fitting, which is either a) no mean parameter, b) constant mean, c) AR(k) structure, or d) general Bayesian regression. In the last case, new regressors need to be provided for prediction.
## S3 method for class 'svdraws' predict(object, steps = 1L, newdata = NULL, ...)## S3 method for class 'svdraws' predict(object, steps = 1L, newdata = NULL, ...)
object |
|
steps |
optional single number, coercible to integer. Denotes the number of steps to forecast. |
newdata |
only in case d) of the description corresponds to input
parameter |
... |
currently ignored. |
Returns an object of class svpredict, a list containing
three elements:
vol |
|
h |
|
y |
|
You can use the resulting object within plot.svdraws (see example below), or use
the list items in the usual coda methods for mcmc objects to
print, plot, or summarize the predictions.
# 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$y, draws = 5000, burnin = 100, priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2) ## Predict 10 days ahead fore <- predict(draws, 10) ## Check out the results summary(predlatent(fore)) summary(predy(fore)) plot(draws, forecast = fore) # Example 2 ## Simulate now an SV process with an AR(1) mean structure len <- 109L simar <- svsim(len, phi = 0.93, sigma = 0.15, mu = -9) for (i in 2:len) { simar$y[i] <- 0.1 - 0.7 * simar$y[i-1] + simar$vol[i] * rnorm(1) } ## Obtain 7000 draws drawsar <- svsample(simar$y, draws = 7000, burnin = 300, designmatrix = "ar1", priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2) ## Predict 7 days ahead (using AR(1) mean for the returns) forear <- predict(drawsar, 7) ## Check out the results plot(forear) plot(drawsar, forecast = forear) ## Not run: # Example 3 ## Simulate now an SV process with leverage and with non-zero mean len <- 96L regressors <- cbind(rep_len(1, len), rgamma(len, 0.5, 0.25)) betas <- rbind(-1.1, 2) simreg <- svsim(len, rho = -0.42) simreg$y <- simreg$y + as.numeric(regressors %*% betas) ## Obtain 12000 draws drawsreg <- svsample(simreg$y, draws = 12000, burnin = 3000, designmatrix = regressors, priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2, priorrho = c(4, 4)) ## Predict 5 days ahead using new regressors predlen <- 5L predregressors <- cbind(rep_len(1, predlen), rgamma(predlen, 0.5, 0.25)) forereg <- predict(drawsreg, predlen, predregressors) ## Check out the results summary(predlatent(forereg)) summary(predy(forereg)) plot(forereg) plot(drawsreg, forecast = forereg) ## End(Not run)# 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$y, draws = 5000, burnin = 100, priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2) ## Predict 10 days ahead fore <- predict(draws, 10) ## Check out the results summary(predlatent(fore)) summary(predy(fore)) plot(draws, forecast = fore) # Example 2 ## Simulate now an SV process with an AR(1) mean structure len <- 109L simar <- svsim(len, phi = 0.93, sigma = 0.15, mu = -9) for (i in 2:len) { simar$y[i] <- 0.1 - 0.7 * simar$y[i-1] + simar$vol[i] * rnorm(1) } ## Obtain 7000 draws drawsar <- svsample(simar$y, draws = 7000, burnin = 300, designmatrix = "ar1", priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2) ## Predict 7 days ahead (using AR(1) mean for the returns) forear <- predict(drawsar, 7) ## Check out the results plot(forear) plot(drawsar, forecast = forear) ## Not run: # Example 3 ## Simulate now an SV process with leverage and with non-zero mean len <- 96L regressors <- cbind(rep_len(1, len), rgamma(len, 0.5, 0.25)) betas <- rbind(-1.1, 2) simreg <- svsim(len, rho = -0.42) simreg$y <- simreg$y + as.numeric(regressors %*% betas) ## Obtain 12000 draws drawsreg <- svsample(simreg$y, draws = 12000, burnin = 3000, designmatrix = regressors, priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2, priorrho = c(4, 4)) ## Predict 5 days ahead using new regressors predlen <- 5L predregressors <- cbind(rep_len(1, predlen), rgamma(predlen, 0.5, 0.25)) forereg <- predict(drawsreg, predlen, predregressors) ## Check out the results summary(predlatent(forereg)) summary(predy(forereg)) plot(forereg) plot(drawsreg, forecast = forereg) ## End(Not run)
This function gives access to a larger set of prior distributions in case the default choice is unsatisfactory.
specify_priors( mu = sv_normal(mean = 0, sd = 100), phi = sv_beta(shape1 = 5, shape2 = 1.5), sigma2 = sv_gamma(shape = 0.5, rate = 0.5), nu = sv_infinity(), rho = sv_constant(0), latent0_variance = "stationary", beta = sv_multinormal(mean = 0, sd = 10000, dim = 1) )specify_priors( mu = sv_normal(mean = 0, sd = 100), phi = sv_beta(shape1 = 5, shape2 = 1.5), sigma2 = sv_gamma(shape = 0.5, rate = 0.5), nu = sv_infinity(), rho = sv_constant(0), latent0_variance = "stationary", beta = sv_multinormal(mean = 0, sd = 10000, dim = 1) )
mu |
one of |
phi |
one of |
sigma2 |
one of |
nu |
one of |
rho |
one of |
latent0_variance |
either the character string |
beta |
an |
Other priors:
sv_constant()
stochvol
The functions below can be supplied to specify_priors
to overwrite the default set of prior distributions in svsample.
The functions have mean, range, density, and
print methods.
sv_constant(value) sv_normal(mean = 0, sd = 1) sv_multinormal(mean = 0, precision = NULL, sd = 1, dim = NA) sv_gamma(shape, rate) sv_inverse_gamma(shape, scale) sv_beta(shape1, shape2) sv_exponential(rate) sv_infinity()sv_constant(value) sv_normal(mean = 0, sd = 1) sv_multinormal(mean = 0, precision = NULL, sd = 1, dim = NA) sv_gamma(shape, rate) sv_inverse_gamma(shape, scale) sv_beta(shape1, shape2) sv_exponential(rate) sv_infinity()
value |
The constant value for the degenerate constant distribution |
mean |
Expected value for the univariate normal distribution or mean vector of the multivariate normal distribution |
sd |
Standard deviation for the univariate normal distribution or constant scale of the multivariate normal distribution |
precision |
Precision matrix for the multivariate normal distribution |
dim |
(optional) Dimension of the multivariate distribution |
shape |
Shape parameter for the distribution |
rate |
Rate parameter for the distribution |
scale |
Scale parameter for the distribution |
shape1 |
First shape parameter for the distribution |
shape2 |
Second shape parameter for the distribution |
Multivariate normal objects can be specified several ways. The most general way is by calling
sv_multinormal(mean, precision), which allows for arbitrary mean and (valid) precision
arguments. Constant mean vectors and constant diagonal precision matrices of dimension D
can be created two ways: either sv_multinormal(mean, sd, dim = D) or
rep(sv_normal(mean, sd), length.out = D).
Other priors:
specify_priors()
svlm is a wrapper around svsample with a formula interface.
The name derives from SV and lm because a linear model with SV residuals is fitted.
The function simulates from the joint posterior distribution of the regression coefficients and 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.
svlm( formula, data, draws = 10000, burnin = 1000, heavytails = FALSE, asymmetry = FALSE, priorspec = NULL, thin = 1, 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, ... )svlm( formula, data, draws = 10000, burnin = 1000, heavytails = FALSE, asymmetry = FALSE, priorspec = NULL, thin = 1, 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, ... )
formula |
an object of class |
data |
an optional data frame, list or environment (or object
coercible by |
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. |
heavytails |
if |
asymmetry |
if |
priorspec |
using the smart constructor |
thin |
single number greater or equal to 1, coercible to integer.
Every |
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
|
startpara |
optional named list, containing the starting values
for the parameter draws. If supplied, |
startlatent |
optional vector of length |
parallel |
optional one of |
n_cpus |
optional positive integer, the number of CPUs to be used in case of
parallel computations. Defaults to |
cl |
optional so-called SNOW cluster object as implemented in package
|
n_chains |
optional positive integer specifying the number of independent MCMC chains |
print_progress |
optional one of |
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
|
... |
Any extra arguments will be forwarded to
|
For details concerning the algorithm please see the paper by Kastner and Frühwirth-Schnatter (2014) and Hosszejni and Kastner (2019).
The value returned is a list object of class svdraws holding
para |
|
latent |
|
latent0 |
|
tau |
|
beta |
|
y |
the left hand side of the observation equation, usually
the argument |
runtime |
|
priors |
a |
thinning |
|
summary |
|
meanmodel |
|
svlm |
a flag for the use of |
model_terms |
helper object that represents the formula |
formula |
argument |
xlevels |
helper object that is needed to interpret the formula |
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.
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.
svsample, svsim, specify_priors
# Simulate data n <- 50L dat <- data.frame(x = runif(n, 3, 4), z = runif(n, -1, -0.5)) designmatrix <- matrix(c(dat$x, dat$x^2, log10(dat$x), dat$z), ncol = 4) betas <- matrix(c(-1, 1, 2, 0), ncol = 1) y <- designmatrix %*% betas + svsim(n)$y dat$y <- y # Formula interface res <- svlm(y ~ 0 + x + I(x^2) + log10(x) + z, data = dat) # Prediction predn <- 10L preddat <- data.frame(x = runif(predn, 3, 4), z = runif(predn, -1, -0.5)) pred <- predict(res, newdata = preddat, steps = predn)# Simulate data n <- 50L dat <- data.frame(x = runif(n, 3, 4), z = runif(n, -1, -0.5)) designmatrix <- matrix(c(dat$x, dat$x^2, log10(dat$x), dat$z), ncol = 4) betas <- matrix(c(-1, 1, 2, 0), ncol = 1) y <- designmatrix %*% betas + svsim(n)$y dat$y <- y # Formula interface res <- svlm(y ~ 0 + x + I(x^2) + log10(x) + z, data = dat) # Prediction predn <- 10L preddat <- data.frame(x = runif(predn, 3, 4), z = runif(predn, -1, -0.5)) pred <- predict(res, newdata = preddat, steps = predn)
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.
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 )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 )
y |
numeric vector containing the data (usually log-returns), which
must not contain zeros. Alternatively, |
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 |
priormu |
numeric vector of length 2, indicating mean and standard
deviation for the Gaussian prior distribution of the parameter |
priorphi |
numeric vector of length 2, indicating the shape parameters
for the Beta prior distribution of the transformed parameter
|
priorsigma |
single positive real number, which stands for the scaling
of the transformed parameter |
priornu |
single non-negative number, indicating the rate parameter
for the exponential prior distribution of the parameter; can be |
priorrho |
either |
priorbeta |
numeric vector of length 2, indicating the mean and
standard deviation of the Gaussian prior for the regression parameters. The
default value is |
priorlatent0 |
either a single non-negative number or the string
|
priorspec |
in case one needs different prior distributions than the
ones specified by |
thin |
single number greater or equal to 1, coercible to integer.
Every |
thinpara |
single number greater or equal to 1, coercible to integer.
Every |
thinlatent |
single number greater or equal to 1, coercible to integer.
Every |
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
|
startpara |
optional named list, containing the starting values
for the parameter draws. If supplied, |
startlatent |
optional vector of length |
parallel |
optional one of |
n_cpus |
optional positive integer, the number of CPUs to be used in case of
parallel computations. Defaults to |
cl |
optional so-called SNOW cluster object as implemented in package
|
n_chains |
optional positive integer specifying the number of independent MCMC chains |
print_progress |
optional one of |
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
|
... |
Any extra arguments will be forwarded to
|
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).
The value returned is a list object of class svdraws holding
para |
|
latent |
|
latent0 |
|
tau |
|
beta |
|
y |
the left hand side of the observation equation, usually
the argument |
runtime |
|
priors |
a |
thinning |
|
summary |
|
meanmodel |
|
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.
If y contains zeros, you might want to consider de-meaning your
returns or use designmatrix = "ar0".
svsample2 is deprecated.
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.
############### # 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)############### # 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)
C++ Functions in stochvol
All the heavy lifting in stochvol is implemented in C++
with the help of R packages Rcpp and RcppArmadillo.
These functions call the MCMC samplers in C++ directly without any
any validation and transformations, expert use only!
svsample_fast_cpp( y, draws = 1, burnin = 0, designmatrix = matrix(NA), priorspec = specify_priors(), thinpara = 1, thinlatent = 1, keeptime = "all", startpara, startlatent, keeptau = !inherits(priorspec$nu, "sv_infinity"), print_settings = list(quiet = TRUE, n_chains = 1, chain = 1), correct_model_misspecification = FALSE, interweave = TRUE, myoffset = 0, fast_sv = get_default_fast_sv() ) svsample_general_cpp( y, draws = 1, burnin = 0, designmatrix = matrix(NA), priorspec = specify_priors(), thinpara = 1, thinlatent = 1, keeptime = "all", startpara, startlatent, keeptau = !inherits(priorspec$nu, "sv_infinity"), print_settings = list(quiet = TRUE, n_chains = 1, chain = 1), correct_model_misspecification = FALSE, interweave = TRUE, myoffset = 0, general_sv = get_default_general_sv(priorspec) )svsample_fast_cpp( y, draws = 1, burnin = 0, designmatrix = matrix(NA), priorspec = specify_priors(), thinpara = 1, thinlatent = 1, keeptime = "all", startpara, startlatent, keeptau = !inherits(priorspec$nu, "sv_infinity"), print_settings = list(quiet = TRUE, n_chains = 1, chain = 1), correct_model_misspecification = FALSE, interweave = TRUE, myoffset = 0, fast_sv = get_default_fast_sv() ) svsample_general_cpp( y, draws = 1, burnin = 0, designmatrix = matrix(NA), priorspec = specify_priors(), thinpara = 1, thinlatent = 1, keeptime = "all", startpara, startlatent, keeptau = !inherits(priorspec$nu, "sv_infinity"), print_settings = list(quiet = TRUE, n_chains = 1, chain = 1), correct_model_misspecification = FALSE, interweave = TRUE, myoffset = 0, general_sv = get_default_general_sv(priorspec) )
y |
numeric vector of the observations |
draws |
single positive integer, the number of draws to return (after the burn-in) |
burnin |
single positive integer, length of warm-up period, this number of draws are discarded from the beginning |
designmatrix |
numeric matrix of covariates. Dimensions:
|
priorspec |
a |
thinpara |
single number greater or equal to 1, coercible to integer.
Every |
thinlatent |
single number greater or equal to 1, coercible to integer.
Every |
keeptime |
Either 'all' (the default) or 'last'. Indicates which latent volatility draws should be stored. |
startpara |
named list, containing the starting values for the parameter draws. It must contain elements
|
startlatent |
vector of length |
keeptau |
Logical value indicating whether the 'variance inflation factors' should be stored (used for the sampler with conditional t innovations only). This may be useful to check at what point(s) in time the normal disturbance had to be 'upscaled' by a mixture factor and when the series behaved 'normally'. |
print_settings |
List of three elements:
Please note that this function does not run multiple independent chains
but |
correct_model_misspecification |
Logical value. If |
interweave |
Logical value. If |
myoffset |
Single non-negative number that is used in
|
fast_sv |
named list of expert settings. We recommend the use of |
general_sv |
named list of expert settings. We recommend the use of |
The sampling functions are separated into fast SV and general SV. See more details in the sections below.
Fast SV was developed in Kastner and Fruehwirth-Schnatter (2014). Fast SV estimates an approximate SV model without leverage, where the approximation comes in through auxiliary mixture approximations to the exact SV model. The sampler uses the ancillarity-sufficiency interweaving strategy (ASIS) to improve on the sampling efficiency of the model parameters, and it employs all-without-a-loop (AWOL) for computationally efficient Kalman filtering of the conditionally Gaussian state space. Correction for model misspecification happens as a post-processing step.
Fast SV employs sampling strategies that have been fine-tuned and specified for
vanilla SV (no leverage), and hence it can be fast and efficient but also more limited
in its feature set. The conditions for the fast SV sampler: rho == 0; mu
has either a normal prior or it is also constant 0; the prior for phi
is a beta distribution; the prior for sigma^2 is either a gamma distribution
with shape 0.5 or a mean- and variance-matched inverse gamma distribution;
either keeptime == 'all' or correct_model_misspecification == FALSE.
These criteria are NOT VALIDATED by fast SV on the C++ level!
General SV also estimates an approximate SV model without leverage, where the approximation comes in through auxiliary mixture approximations to the exact SV model. The sampler uses both ASIS and AWOL.
General SV employs adapted random walk Metropolis-Hastings as the proposal for
the parameters mu, phi, sigma, and rho. Therefore,
more general prior distributions are allowed in this case.
# Draw one sample using fast SV and general SV y <- svsim(40)$y params <- list(mu = -10, phi = 0.9, sigma = 0.1, nu = Inf, rho = 0, beta = NA, latent0 = -10) res_fast <- svsample_fast_cpp(y, startpara = params, startlatent = rep(-10, 40)) res_gen <- svsample_general_cpp(y, startpara = params, startlatent = rep(-10, 40)) # Embed SV in another sampling scheme ## vanilla SV len <- 40L draws <- 1000L burnin <- 200L param_store <- matrix(NA, draws, 3, dimnames = list(NULL, c("mu", "phi", "sigma"))) startpara <- list(mu = 0, phi = 0.9, sigma = 0.1, nu = Inf, rho = 0, beta = NA, latent0 = 0) startlatent <- rep(0, len) for (i in seq_len(burnin+draws)) { # draw the data in the bigger sampling scheme # now we simulate y from vanilla SV y <- svsim(len, mu = 0, phi = 0.9, sigma = 0.1)$y # call SV sampler res <- svsample_fast_cpp(y, startpara = startpara, startlatent = startlatent) # administrate values startpara[c("mu","phi","sigma")] <- as.list(res$para[, c("mu", "phi", "sigma")]) startlatent <- drop(res$latent) # store draws after the burnin if (i > burnin) { param_store[i-burnin, ] <- res$para[, c("mu", "phi", "sigma")] } } ### quick look at the traceplots ts.plot(param_store, col = 1:3)# Draw one sample using fast SV and general SV y <- svsim(40)$y params <- list(mu = -10, phi = 0.9, sigma = 0.1, nu = Inf, rho = 0, beta = NA, latent0 = -10) res_fast <- svsample_fast_cpp(y, startpara = params, startlatent = rep(-10, 40)) res_gen <- svsample_general_cpp(y, startpara = params, startlatent = rep(-10, 40)) # Embed SV in another sampling scheme ## vanilla SV len <- 40L draws <- 1000L burnin <- 200L param_store <- matrix(NA, draws, 3, dimnames = list(NULL, c("mu", "phi", "sigma"))) startpara <- list(mu = 0, phi = 0.9, sigma = 0.1, nu = Inf, rho = 0, beta = NA, latent0 = 0) startlatent <- rep(0, len) for (i in seq_len(burnin+draws)) { # draw the data in the bigger sampling scheme # now we simulate y from vanilla SV y <- svsim(len, mu = 0, phi = 0.9, sigma = 0.1)$y # call SV sampler res <- svsample_fast_cpp(y, startpara = startpara, startlatent = startlatent) # administrate values startpara[c("mu","phi","sigma")] <- as.list(res$para[, c("mu", "phi", "sigma")]) startlatent <- drop(res$latent) # store draws after the burnin if (i > burnin) { param_store[i-burnin, ] <- res$para[, c("mu", "phi", "sigma")] } } ### quick look at the traceplots ts.plot(param_store, col = 1:3)
svsample_roll performs rolling window estimation based on svsample.
A convenience function for backtesting purposes.
svsample_roll( y, designmatrix = NA, n_ahead = 1, forecast_length = 500, n_start = NULL, refit_every = 1, refit_window = c("moving", "expanding"), calculate_quantile = c(0.01), calculate_predictive_likelihood = TRUE, keep_draws = FALSE, parallel = c("no", "multicore", "snow"), n_cpus = 1L, cl = NULL, ... ) svtsample_roll( y, designmatrix = NA, n_ahead = 1, forecast_length = 500, n_start = NULL, refit_every = 1, refit_window = c("moving", "expanding"), calculate_quantile = c(0.01), calculate_predictive_likelihood = TRUE, keep_draws = FALSE, parallel = c("no", "multicore", "snow"), n_cpus = 1L, cl = NULL, ... ) svlsample_roll( y, designmatrix = NA, n_ahead = 1, forecast_length = 500, n_start = NULL, refit_every = 1, refit_window = c("moving", "expanding"), calculate_quantile = c(0.01), calculate_predictive_likelihood = TRUE, keep_draws = FALSE, parallel = c("no", "multicore", "snow"), n_cpus = 1L, cl = NULL, ... ) svtlsample_roll( y, designmatrix = NA, n_ahead = 1, forecast_length = 500, n_start = NULL, refit_every = 1, refit_window = c("moving", "expanding"), calculate_quantile = c(0.01), calculate_predictive_likelihood = TRUE, keep_draws = FALSE, parallel = c("no", "multicore", "snow"), n_cpus = 1L, cl = NULL, ... )svsample_roll( y, designmatrix = NA, n_ahead = 1, forecast_length = 500, n_start = NULL, refit_every = 1, refit_window = c("moving", "expanding"), calculate_quantile = c(0.01), calculate_predictive_likelihood = TRUE, keep_draws = FALSE, parallel = c("no", "multicore", "snow"), n_cpus = 1L, cl = NULL, ... ) svtsample_roll( y, designmatrix = NA, n_ahead = 1, forecast_length = 500, n_start = NULL, refit_every = 1, refit_window = c("moving", "expanding"), calculate_quantile = c(0.01), calculate_predictive_likelihood = TRUE, keep_draws = FALSE, parallel = c("no", "multicore", "snow"), n_cpus = 1L, cl = NULL, ... ) svlsample_roll( y, designmatrix = NA, n_ahead = 1, forecast_length = 500, n_start = NULL, refit_every = 1, refit_window = c("moving", "expanding"), calculate_quantile = c(0.01), calculate_predictive_likelihood = TRUE, keep_draws = FALSE, parallel = c("no", "multicore", "snow"), n_cpus = 1L, cl = NULL, ... ) svtlsample_roll( y, designmatrix = NA, n_ahead = 1, forecast_length = 500, n_start = NULL, refit_every = 1, refit_window = c("moving", "expanding"), calculate_quantile = c(0.01), calculate_predictive_likelihood = TRUE, keep_draws = FALSE, parallel = c("no", "multicore", "snow"), n_cpus = 1L, cl = NULL, ... )
y |
numeric vector containing the data (usually log-returns), which
must not contain zeros. Alternatively, |
designmatrix |
regression design matrix for modeling the mean. Must
have |
n_ahead |
number of time steps to predict from each time window. |
forecast_length |
the time horizon at the end of the data set that is used for backtesting. |
n_start |
optional the starting time point for backtesting.
Computed from |
refit_every |
the SV model is refit every |
refit_window |
one of |
calculate_quantile |
vector of numbers between 0 and 1.
These quantiles are predicted using |
calculate_predictive_likelihood |
boolean. If |
keep_draws |
boolean. If |
parallel |
one of |
n_cpus |
optional positive integer, the number of CPUs to be used in case of
parallel computations. Defaults to |
cl |
optional so-called SNOW cluster object as implemented in package
|
... |
Any extra arguments will be forwarded to
|
Functions svtsample_roll, svlsample_roll, and svtlsample_roll are
wrappers around svsample_roll with convenient default values for the SV
model with t-errors, leverage, and both t-errors and leverage, respectively.
The value returned is a list object of class svdraws_roll
holding a list item for every time window. The elements of these list items are
indices |
a list object containing two elements: |
quantiles |
the input parameter |
refit_every |
the input parameter |
predictive_likelihood |
present only if |
predictive_quantile |
present only if |
fit |
present only if |
prediction |
present only if |
To display the output, use print and summary. The
print method simply prints a short summary of the setup;
the summary method displays the summary statistics
of the backtesting.
The function executes svsample (length(y) - arorder - n_ahead - n_start + 2) %/% refit_every times.
svsim, specify_priors, svsample
# Simulate from the true model sim <- svsim(200) # Perform rolling estimation using the vanilla SV # model and default priors roll <- svsample_roll(sim, draws = 5000, burnin = 2000, keep_draws = TRUE, forecast_length = 10, n_ahead = 1, refit_every = 1, refit_window = "moving", calculate_predictive_likelihood = TRUE, calculate_quantile = c(0.01, 0.05)) # Perform rolling estimation by making use # of two CPU cores, advanced priors, and multiple # chains with pre-set initial values. Let us combine # that with an AR(2) specification prior_beta <- sv_multinormal(c(1,0,-1), rbind(c(1, 0, 0.1), c(0, 0.3, -0.04), c(0.1, -0.04, 0.1))) priorspec <- specify_priors(rho = sv_beta(4, 4), latent0_variance = sv_constant(1), beta = prior_beta, nu = sv_exponential(0.05)) startpara <- list(list(mu = -9, phi = 0.3), list(mu = -11, sigma = 0.1, phi = 0.95), list(phi = 0.99)) roll <- svsample_roll(sim, draws = 5000, burnin = 2000, designmatrix = "ar2", priorspec = priorspec, startpara = startpara, parallel = "snow", n_cpus = 2, n_chains = 3, keep_draws = TRUE, forecast_length = 10, n_ahead = 1, refit_every = 1, refit_window = "expanding", calculate_predictive_likelihood = TRUE, calculate_quantile = c(0.01, 0.05))# Simulate from the true model sim <- svsim(200) # Perform rolling estimation using the vanilla SV # model and default priors roll <- svsample_roll(sim, draws = 5000, burnin = 2000, keep_draws = TRUE, forecast_length = 10, n_ahead = 1, refit_every = 1, refit_window = "moving", calculate_predictive_likelihood = TRUE, calculate_quantile = c(0.01, 0.05)) # Perform rolling estimation by making use # of two CPU cores, advanced priors, and multiple # chains with pre-set initial values. Let us combine # that with an AR(2) specification prior_beta <- sv_multinormal(c(1,0,-1), rbind(c(1, 0, 0.1), c(0, 0.3, -0.04), c(0.1, -0.04, 0.1))) priorspec <- specify_priors(rho = sv_beta(4, 4), latent0_variance = sv_constant(1), beta = prior_beta, nu = sv_exponential(0.05)) startpara <- list(list(mu = -9, phi = 0.3), list(mu = -11, sigma = 0.1, phi = 0.95), list(phi = 0.99)) roll <- svsample_roll(sim, draws = 5000, burnin = 2000, designmatrix = "ar2", priorspec = priorspec, startpara = startpara, parallel = "snow", n_cpus = 2, n_chains = 3, keep_draws = TRUE, forecast_length = 10, n_ahead = 1, refit_every = 1, refit_window = "expanding", calculate_predictive_likelihood = TRUE, calculate_quantile = c(0.01, 0.05))
svsim is used to produce realizations of a stochastic volatility (SV)
process.
svsim(len, mu = -10, phi = 0.98, sigma = 0.2, nu = Inf, rho = 0)svsim(len, mu = -10, phi = 0.98, sigma = 0.2, nu = Inf, rho = 0)
len |
length of the simulated time series. |
mu |
level of the latent log-volatility AR(1) process. The defaults
value is |
phi |
persistence of the latent log-volatility AR(1) process. The
default value is |
sigma |
volatility of the latent log-volatility AR(1) process. The
default value is |
nu |
degrees-of-freedom for the conditional innovations distribution.
The default value is |
rho |
correlation between the observation and the increment of the
log-volatility. The default value is |
This function draws an initial log-volatility h_0 from the stationary
distribution of the AR(1) process defined by phi, sigma, and mu.
Then the function jointly simulates the log-volatility series
h_1,...,h_n with the given AR(1) structure, and the “log-return” series
y_1,...,y_n with mean 0 and standard deviation exp(h/2).
Additionally, for each index i, y_i can be set to have a conditionally heavy-tailed
residual (through nu) and/or to be correlated with (h_{i+1}-h_i)
(through rho, the so-called leverage effect, resulting in asymmetric “log-returns”).
The output is a list object of class svsim containing
vector of length len containing the simulated data,
usually interpreted as “log-returns”.
vector of length
len containing the simulated instantaneous volatilities.
These are if nu == Inf, and they are
for finite nu.
The initial volatility exp(h_0/2),
drawn from the stationary distribution of the latent AR(1) process.
a named list with five elements mu, phi,
sigma, nu, and rho, containing
the corresponding arguments.
vector of the latent state space for .
initial element of the latent state space .
vector of length len containing the simulated auxiliary
variables for the Student-t residuals when nu is finite. More precisely,
.
The function generates the “log-returns” by
y <- exp(-h/2)*rt(h, df=nu). That means that in the case of nu < Inf
the (conditional) volatility is sqrt(nu/(nu-2))*exp(h/2), and that corrected value
is shown in the print, summary and plot methods.
To display the output use print, summary and plot. The
print method simply prints the content of the object in a moderately
formatted manner. The summary method provides some summary statistics
(in %), and the plot method plots the the simulated 'log-returns'
y along with the corresponding volatilities vol.
Gregor Kastner [email protected]
## Simulate a highly persistent SV process of length 500 sim <- svsim(500, phi = 0.99, sigma = 0.1) print(sim) summary(sim) plot(sim) ## Simulate an SV process with leverage sim <- svsim(200, phi = 0.94, sigma = 0.15, rho = -0.6) print(sim) summary(sim) plot(sim) ## Simulate an SV process with conditionally heavy-tails sim <- svsim(250, phi = 0.91, sigma = 0.05, nu = 5) print(sim) summary(sim) plot(sim)## Simulate a highly persistent SV process of length 500 sim <- svsim(500, phi = 0.99, sigma = 0.1) print(sim) summary(sim) plot(sim) ## Simulate an SV process with leverage sim <- svsim(200, phi = 0.94, sigma = 0.15, rho = -0.6) print(sim) summary(sim) plot(sim) ## Simulate an SV process with conditionally heavy-tails sim <- svsim(250, phi = 0.91, sigma = 0.05, nu = 5) print(sim) summary(sim) plot(sim)
Samples the mixture indicators, the latent variables, and the model independent parameters mu, phi, and sigma. The input is the logarithm of the squared de-meaned observations. An approximate SV model is estimated instead of the exact SV model by the use of auxiliary mixture sampling. Depending on the prior specification, mu might not be updated. Depending on the expert settings, the function might follow the ancillarity-sufficiency interweaving strategy (ASIS, Yu and Meng, 2011) for sampling mu, phi, and sigma. Furthermore, the user can turn off the sampling of the parameters, the latents, or the mixture indicators in the expert settings.
update_fast_sv(log_data2, mu, phi, sigma, h0, h, r, prior_spec, expert)update_fast_sv(log_data2, mu, phi, sigma, h0, h, r, prior_spec, expert)
log_data2 |
log(data^2), where data is the vector of de-meaned observations |
mu |
parameter mu. Level of the latent process h. Updated in place |
phi |
parameter phi, persistence of the latent process h. Updated in place |
sigma |
parameter sigma, volatility of the latent process h, also called volvol. Updated in place |
h0 |
parameter h0, the initial value of the latent process h. Updated in place |
h |
the vector of the latent process. Updated in place |
r |
the vector of the mixture indicators. Updated in place |
prior_spec |
prior specification object. See type_definitions.h |
expert |
expert settings for this function. See type_definitions.h |
Other stochvol_cpp:
update_general_sv(),
update_regressors(),
update_t_error()
Samples the latent variables and the model independent parameters mu, phi, sigma, and rho. The observations need to be provided in different formats for efficiency. An approximate SV model is as the default posterior distribution for the latent vector; however, there is the option to correct for model misspecification through the expert settings. Depending on the prior specification, some of mu, phi, sigma, and rho might not be updated. Depending on the expert settings, the function might follow the ancillarity-sufficiency interweaving strategy (ASIS, Yu and Meng, 2011) for sampling mu, phi, sigma, and rho. Also controlled by the expert settings, Furthermore, the user can turn off the sampling of the parameters, the latents, or the mixture indicators in the expert settings.
update_general_sv( data, log_data2, sign_data, mu, phi, sigma, rho, h0, h, adaptation, prior_spec, expert )update_general_sv( data, log_data2, sign_data, mu, phi, sigma, rho, h0, h, adaptation, prior_spec, expert )
data |
the vector of de-meaned observations |
log_data2 |
log(data^2), where data is the vector of de-meaned observations |
sign_data |
the sign of the data |
mu |
parameter mu. Level of the latent process h. Updated in place |
phi |
parameter phi, persistence of the latent process h. Updated in place |
sigma |
parameter sigma, volatility of the latent process h, also called volvol. Updated in place |
rho |
parameter rho. Accounts for asymmetry/the leverage effect. Updated in place |
h0 |
parameter h0, the initial value of the latent process h. Updated in place |
h |
the vector of the latent process. Updated in place |
adaptation |
object implementing the adaptive Metropolis-Hastings scheme. Updated in place. See adaptation.hpp |
prior_spec |
prior specification object. See type_definitions.h |
expert |
expert settings for this function. See type_definitions.h |
Other stochvol_cpp:
update_fast_sv(),
update_regressors(),
update_t_error()
Samples the coefficients of a linear regression. The prior distribution is multivariate normal and it is specified in prior_spec.
update_regressors(dependent_variable, independent_variables, beta, prior_spec)update_regressors(dependent_variable, independent_variables, beta, prior_spec)
dependent_variable |
the left hand side |
independent_variables |
the matrix of the independent variables. Has to be of same height as the length of the dependent variable |
beta |
the vector of the latent states used in MDA. Updated in place |
prior_spec |
prior specification object. See type_definitions.h |
Other stochvol_cpp:
update_fast_sv(),
update_general_sv(),
update_t_error()
Samples the degrees of freedom parameter of standardized and homoskedastic t-distributed input variates. Marginal data augmentation (MDA) is applied, tau is the vector of auxiliary latent states. Depending on the prior specification, nu might not be updated, just tau.
update_t_error( homosked_data, tau, mean, sd, nu, prior_spec, do_tau_acceptance_rejection = TRUE )update_t_error( homosked_data, tau, mean, sd, nu, prior_spec, do_tau_acceptance_rejection = TRUE )
homosked_data |
de-meaned and homoskedastic observations |
tau |
the vector of the latent states used in MDA. Updated in place |
mean |
the vector of the conditional means // TODO update docs in R |
sd |
the vector of the conditional standard deviations |
nu |
parameter nu. The degrees of freedom for the t-distribution. Updated in place |
prior_spec |
prior specification object. See type_definitions.h |
do_tau_acceptance_rejection |
boolean. If |
The function samples tau and nu from the following hierarchical model: homosked_data_i = sqrt(tau_i) * (mean_i + sd_i * N(0, 1)) tau_i ~ InvGamma(.5*nu, .5*(nu-2)) Naming: The data is homoskedastic ex ante in the model, mean_i and sd_i are conditional on some other parameter in the model. The prior on tau corresponds to a standardized t-distributed heavy tail on the data.
Other stochvol_cpp:
update_fast_sv(),
update_general_sv(),
update_regressors()
Creates or updates a summary of an svdraws object.
updatesummary( x, quantiles = c(0.05, 0.5, 0.95), esspara = TRUE, esslatent = FALSE )updatesummary( x, quantiles = c(0.05, 0.5, 0.95), esspara = TRUE, esslatent = FALSE )
x |
|
quantiles |
numeric vector of posterior quantiles to be computed. The
default is |
esspara |
logical value which indicates whether the effective sample
size (ESS) should be calculated for the parameter draws. This is
achieved by calling |
esslatent |
logical value which indicates whether the effective sample
size (ESS) should be calculated for the latent log-volatility draws.
This is achieved by calling |
updatesummary will always calculate the posterior mean and the
posterior standard deviation of the raw draws and some common
transformations thereof. Moroever, the posterior quantiles, specified by the
argument quantiles, are computed. If esspara and/or
esslatent are TRUE, the corresponding effective sample size
(ESS) will also be included.
The value returned is an updated list object of class svdraws
holding
para |
|
latent |
|
latent0 |
|
y |
argument |
runtime |
|
priors |
|
thinning |
|
summary |
|
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 gives a graphical
overview of the posterior distribution by calling volplot,
traceplot and densplot and displaying the
results on a single page.
updatesummary does not actually overwrite the object's current
summary, but in fact creates a new object with an updated summary. Thus,
don't forget to overwrite the old object if this is want you intend to do.
See the examples below for more details.
## Here is a baby-example to illustrate the idea. ## Simulate an SV time series of length 51 with default parameters: sim <- svsim(51) ## Draw from the posterior: res <- svsample(sim$y, draws = 2000, priorphi = c(10, 1.5)) ## Check out the results: summary(res) plot(res) ## Look at other quantiles and calculate ESS of latents: newquants <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99) res <- updatesummary(res, quantiles = newquants, esslatent = TRUE) ## See the difference? summary(res) plot(res)## Here is a baby-example to illustrate the idea. ## Simulate an SV time series of length 51 with default parameters: sim <- svsim(51) ## Draw from the posterior: res <- svsample(sim$y, draws = 2000, priorphi = c(10, 1.5)) ## Check out the results: summary(res) plot(res) ## Look at other quantiles and calculate ESS of latents: newquants <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99) res <- updatesummary(res, quantiles = newquants, esslatent = TRUE) ## See the difference? summary(res) plot(res)
A helper function that validates the input and extends it with default values if there are missing parts for argument 'expert'.
validate_and_process_expert(expert = NULL, priorspec = specify_priors())validate_and_process_expert(expert = NULL, priorspec = specify_priors())
expert |
list, the input values for expert. |
priorspec |
a |
A list that is the input extended by default values. If the input is invalid, an error is thrown.
Displays quantiles of the posterior distribution of the volatilities over time as well as predictive distributions of future volatilities.
volplot( x, forecast = 0, dates = NULL, show0 = FALSE, forecastlty = NULL, tcl = -0.4, mar = c(1.9, 1.9, 1.9, 0.5), mgp = c(2, 0.6, 0), simobj = NULL, newdata = NULL, ... )volplot( x, forecast = 0, dates = NULL, show0 = FALSE, forecastlty = NULL, tcl = -0.4, mar = c(1.9, 1.9, 1.9, 0.5), mgp = c(2, 0.6, 0), simobj = NULL, newdata = NULL, ... )
x |
|
forecast |
nonnegative integer or object of class |
dates |
vector of length |
show0 |
logical value, indicating whether the initial volatility
|
forecastlty |
vector of line type values (see
|
tcl |
The length of tick marks as a fraction of the height of a line of
text. See |
mar |
numerical vector of length 4, indicating the plot margins. See
|
mgp |
numerical vector of length 3, indicating the axis and label
positions. See |
simobj |
object of class |
newdata |
corresponds to parameter |
... |
further arguments are passed on to the invoked |
Called for its side effects. Returns argument x invisibly.
In case you want different quantiles to be plotted, use
updatesummary on the svdraws object first. An example
of doing so is given below.
Gregor Kastner [email protected]
updatesummary, predict.svdraws
Other plotting:
paradensplot(),
paratraceplot(),
paratraceplot.svdraws(),
plot.svdraws(),
plot.svpredict()
## 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$y, draws = 5000, burnin = 100, priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2) ## Plot the latent volatilities and some forecasts volplot(draws, forecast = 10) ## Re-plot with different quantiles newquants <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99) draws <- updatesummary(draws, quantiles = newquants) volplot(draws, forecast = 10)## 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$y, draws = 5000, burnin = 100, priormu = c(-10, 1), priorphi = c(20, 1.5), priorsigma = 0.2) ## Plot the latent volatilities and some forecasts volplot(draws, forecast = 10) ## Re-plot with different quantiles newquants <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99) draws <- updatesummary(draws, quantiles = newquants) volplot(draws, forecast = 10)