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

Description

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

Usage

svsample(y, draws = 10000, burnin = 1000, designmatrix = NA,
  priormu = c(0, 100), priorphi = c(5, 1.5), priorsigma = 1,
  priornu = NA, priorbeta = c(0, 10000), priorlatent0 = "stationary",
  thinpara = 1, thinlatent = 1, keeptime = "all", thintime = NULL,
  keeptau = FALSE, quiet = FALSE, startpara, startlatent, expert, ...)
svtsample(y, draws = 10000, burnin = 1000, designmatrix = NA,
  priormu = c(0, 100), priorphi = c(5, 1.5), priorsigma = 1,
  priornu = c(2, 50), priorbeta = c(0, 10000),
  priorlatent0 = "stationary", thinpara = 1, thinlatent = 1,
  keeptime = "all", thintime = NULL, keeptau = FALSE,
  quiet = FALSE, startpara, startlatent, expert, ...)

Arguments

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 warning is thrown.

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

numeric vector of length 2 (or NA), indicating the lower and upper bounds for the uniform prior distribution of the parameter nu, the degrees-of-freedom parameter of the conditional innovations t-distribution. The default value is NA, fixing the degrees-of-freedom to infinity. This corresponds to conditional standard normal innovations, the pre-1.1.0 behavior of stochvol.

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.

thinpara

single number greater or equal to 1, coercible to integer. Every thinparath parameter draw is kept and returned. The default value is 1, corresponding to no thinning of the parameter draws i.e. every draw is stored.

thinlatent

single number greater or equal to 1, coercible to integer. Every thinlatentth latent variable draw is kept and returned. The default value is 1, corresponding to no thinning of the latent variable draws, i.e. every draw is kept.

keeptime

Either 'all' (the default) or 'last'. Indicates which latent

thintime

Deprecated. Use 'keeptime' instead.

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'.

quiet

logical value indicating whether the progress bar and other informative output during sampling should be omitted. The default value is FALSE, implying verbose output.

startpara

optional named list, containing the starting values for the parameter draws. If supplied, startpara must contain three elements named mu, phi, and sigma, where mu is an arbitrary numerical value, phi is a real number between -1 and 1, and sigma is a positive real number. Moreover, if priornu is not NA, startpara must also contain an element named nu (the degrees of freedom parameter for the t-innovations). The default value is equal to the prior mean.

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)).

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:

parameterization: Character string equal to "centered", "noncentered", "GIS_C", or "GIS_NC". Defaults to "GIS_C".

mhcontrol: Single numeric value controlling the proposal density of a Metropolis-Hastings (MH) update step when sampling sigma. If mhcontrol is smaller than 0, an independence proposal will be used, while values greater than zero control the stepsize of a log-random-walk proposal. Defaults to -1.

gammaprior: Single logical value indicating whether a Gamma prior for sigma^2 should be used. If set to FALSE, an Inverse Gamma prior is employed. Defaults to TRUE.

truncnormal: Single logical value indicating whether a truncated Gaussian distribution should be used as proposal for draws of phi. If set to FALSE, a regular Gaussian prior is employed and the draw is immediately discarded when values outside the unit ball happen to be drawn. Defaults to FALSE.

mhsteps: Either 1, 2, or 3. Indicates the number of blocks used for drawing from the posterior of the parameters. Defaults to 2.

proposalvar4sigmaphi: Single positive number indicating the conditional prior variance of sigma*phi in the ridge proposal density for sampling (mu, phi). Defaults to 10^8.

proposalvar4sigmatheta: Single positive number indicating the conditional prior variance of sigma*theta in the ridge proposal density for sampling (mu, phi). Defaults to 10^12.

…

Any extra arguments will be forwarded to updatesummary, controlling the type of statistics calculated for the posterior draws.

Value

The value returned is a list object of class svdraws holding

para

mcmc object containing the parameter draws from the posterior distribution.

latent

mcmc object containing the latent instantaneous log-volatility draws from the posterior distribution.

latent0

mcmc object containing the latent initial log-volatility draws from the posterior distribution.

tau

mcmc object containing the latent variance inflation factors for the sampler with conditional t-innovations (optional).

beta

mcmc object containing the regression coefficient draws from the posterior distribution (optional).

the argument y.

runtime

proc_time object containing the run time of the sampler.

priors

list containing the parameter values of the prior distribution, i.e. the arguments priormu, priorphi, priorsigma, and potentially priornu and priorbeta.

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.

Details

For details concerning the algorithm please see the paper by Kastner and Fr<U+00FC>hwirth-Schnatter (2014).

References

Kastner, G. and Fr<U+00FC>hwirth-Schnatter, S. (2014). Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models. Computational Statistics & Data Analysis, 76, 408--423, http://dx.doi.org/10.1016/j.csda.2013.01.002.

Examples

Run this code

# 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)

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


# }
# NOT RUN {
# Example 2
## AR(1) structure for the mean
data(exrates)
len <- 3000
ahead <- 100
y <- head(exrates$USD, len)

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

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

## Calculate predictive quantiles
predquants <- apply(preddraws$y, 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 3
## Predicting USD based on JPY and GBP in the mean
data(exrates)
len <- 3000
ahead <- 30
## Calculate log-returns
logreturns <- apply(exrates[, c("USD", "JPY", "GBP")], 2,
                    function (x) diff(log(x)))
logretUSD <- logreturns[2:(len+1), "USD"]
regressors <- cbind(1, as.matrix(logreturns[1:len, ]))  # lagged by 1 day

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

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

Run the code above in your browser using DataLab

Data Engineering and BI courses are free this week!