# R/wrappers.R In gregorkastner/stochvol: Efficient Bayesian Inference for Stochastic Volatility (SV) Models

#### Documented in svlsample2svsamplesvsample2svtsample

# R wrapper functions for the main MCMC loop

#' Markov Chain Monte Carlo (MCMC) Sampling for the Stochastic Volatility (SV)
#' Model
#'
#' \code{svsample} simulates from the joint posterior distribution of the SV
#' parameters \code{mu}, \code{phi}, \code{sigma} (and potentially \code{nu}),
#' along with the latent log-volatilities \code{h_0,...,h_n} and returns the
#' MCMC draws. If a design matrix is provided, simple Bayesian regression can
#' also be conducted.
#'
#' For details concerning the algorithm please see the paper by Kastner and
#' Frühwirth-Schnatter (2014).
#'
#' @param y numeric vector containing the data (usually log-returns), which
#' must not contain zeros. Alternatively, \code{y} can be an \code{svsim}
#' object. In this case, the returns will be extracted and a warning is thrown.
#' @param 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.
#' @param 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.
#' @param designmatrix regression design matrix for modeling the mean. Must
#' have \code{length(y)} rows. Alternatively, \code{designmatrix} may be a
#' string of the form \code{"arX"}, where \code{X} is a nonnegative integer. To
#' fit a constant mean model, use \code{designmatrix = "ar0"} (which is
#' equivalent to \code{designmatrix = matrix(1, nrow = length(y))}). To fit an
#' AR(1) model, use \code{designmatrix = "ar1"}, and so on. If some elements of
#' \code{designmatrix} are \code{NA}, the mean is fixed to zero (pre-1.2.0
#' behavior of \pkg{stochvol}).
#' @param priormu numeric vector of length 2, indicating mean and standard
#' deviation for the Gaussian prior distribution of the parameter \code{mu},
#' the level of the log-volatility. The default value is \code{c(0, 100)},
#' which constitutes a practically uninformative prior for common exchange rate
#' datasets, stock returns and the like.
#' @param priorphi numeric vector of length 2, indicating the shape parameters
#' for the Beta prior distribution of the transformed parameter
#' \code{(phi + 1) / 2}, where \code{phi} denotes the persistence of the
#' log-volatility. The default value is \code{c(5, 1.5)}, which constitutes a
#' prior that puts some belief in a persistent log-volatility but also
#' encompasses the region where \code{phi} is around 0.
#' @param priorsigma single positive real number, which stands for the scaling
#' of the transformed parameter \code{sigma^2}, where \code{sigma} denotes the
#' volatility of log-volatility. More precisely, \code{sigma^2 ~ priorsigma *
#' chisq(df = 1)}. The default value is \code{1}, which constitutes a
#' reasonably vague prior for many common exchange rate datasets, stock returns
#' and the like.
#' @param priornu numeric vector of length 2 (or \code{NA}), indicating the
#' lower and upper bounds for the uniform prior distribution of the parameter
#' \code{nu}, the degrees-of-freedom parameter of the conditional innovations
#' t-distribution. The default value is \code{NA}, fixing the
#' degrees-of-freedom to infinity. This corresponds to conditional standard
#' normal innovations, the pre-1.1.0 behavior of \pkg{stochvol}.
#' @param priorbeta numeric vector of length 2, indicating the mean and
#' standard deviation of the Gaussian prior for the regression parameters. The
#' default value is \code{c(0, 10000)}, which constitutes a very vague prior
#' for many common datasets. Not used if \code{designmatrix} is \code{NA}.
#' @param priorlatent0 either a single non-negative number or the string
#' \code{'stationary'} (the default, also the behavior before version 1.3.0).
#' When \code{priorlatent0} is equal to \code{'stationary'}, the stationary
#' distribution of the latent AR(1)-process is used as the prior for the
#' initial log-volatility \code{h_0}. When \code{priorlatent0} is equal to a
#' number \eqn{B}, we have \eqn{h_0 \sim N(\mu, B\sigma^2)} a priori.
#' @param thinpara single number greater or equal to 1, coercible to integer.
#' Every \code{thinpara}th 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.
#' @param thinlatent single number greater or equal to 1, coercible to integer.
#' Every \code{thinlatent}th 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.
#' @param thintime \emph{Deprecated.} Use 'keeptime' instead.
#' @param keeptime Either 'all' (the default) or 'last'. Indicates which latent
#  volatility draws should be stored.
#' @param 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'.
#' @param quiet logical value indicating whether the progress bar and other
#' informative output during sampling should be omitted. The default value is
#' \code{FALSE}, implying verbose output.
#' @param startpara \emph{optional} named list, containing the starting values
#' for the parameter draws. If supplied, \code{startpara} must contain three
#' elements named \code{mu}, \code{phi}, and \code{sigma}, where \code{mu} is
#' an arbitrary numerical value, \code{phi} is a real number between \code{-1}
#' and \code{1}, and \code{sigma} is a positive real number. Moreover, if
#' \code{priornu} is not \code{NA}, \code{startpara} must also contain an
#' element named \code{nu} (the degrees of freedom parameter for the
#' t-innovations). The default value is equal to the prior mean.
#' @param startlatent \emph{optional} vector of length \code{length(y)},
#' containing the starting values for the latent log-volatility draws. The
#' default value is \code{rep(-10, length(y))}.
#' @param expert \emph{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
#' \code{expert} is provided, it may contain the following named elements:
#'
#' \code{parameterization}: Character string equal to \code{"centered"},
#' \code{"noncentered"}, \code{"GIS_C"}, or \code{"GIS_NC"}. Defaults to
#' \code{"GIS_C"}.
#'
#' \code{mhcontrol}: Single numeric value controlling the proposal density of a
#' Metropolis-Hastings (MH) update step when sampling \code{sigma}. If
#' \code{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 \code{-1}.
#'
#' \code{gammaprior}: Single logical value indicating whether a Gamma prior for
#' \code{sigma^2} should be used. If set to \code{FALSE}, an Inverse Gamma
#' prior is employed. Defaults to \code{TRUE}.
#'
#' \code{truncnormal}: Single logical value indicating whether a truncated
#' Gaussian distribution should be used as proposal for draws of \code{phi}. If
#' set to \code{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 \code{FALSE}.
#'
#' \code{mhsteps}: Either \code{1}, \code{2}, or \code{3}. Indicates the number
#' of blocks used for drawing from the posterior of the parameters. Defaults to
#' \code{2}.
#'
#' \code{proposalvar4sigmaphi}: Single positive number indicating the
#' conditional prior variance of \code{sigma*phi} in the ridge \emph{proposal}
#' density for sampling \code{(mu, phi)}. Defaults to \code{10^8}.
#'
#' \code{proposalvar4sigmatheta}: Single positive number indicating the
#' conditional prior variance of \code{sigma*theta} in the ridge
#' \emph{proposal} density for sampling \code{(mu, phi)}. Defaults to
#' \code{10^12}.
#' @param \dots Any extra arguments will be forwarded to
#' \code{\link{updatesummary}}, controlling the type of statistics calculated
#' for the posterior draws.
#' @return The value returned is a list object of class \code{svdraws} holding
#' \item{para}{\code{mcmc} object containing the \emph{parameter} draws from
#' the posterior distribution.}
#' \item{latent}{\code{mcmc} object containing the
#' \emph{latent instantaneous log-volatility} draws from the posterior
#' distribution.}
#' \item{latent0}{\code{mcmc} object containing the \emph{latent
#' initial log-volatility} draws from the posterior distribution.}
#' \item{tau}{\code{mcmc} object containing the \emph{latent variance inflation
#' factors} for the sampler with conditional t-innovations \emph{(optional)}.}
#' \item{beta}{\code{mcmc} object containing the \emph{regression coefficient}
#' draws from the posterior distribution \emph{(optional)}.}
#' \item{y}{the
#' argument \code{y}.}
#' \item{runtime}{\code{proc_time} object containing the
#' run time of the sampler.}
#' \item{priors}{\code{list} containing the parameter
#' values of the prior distribution, i.e. the arguments \code{priormu},
#' \code{priorphi}, \code{priorsigma}, and potentially \code{priornu} and
#' \code{priorbeta}.}
#' \item{thinning}{\code{list} containing the thinning
#' parameters, i.e. the arguments \code{thinpara}, \code{thinlatent} and
#' \code{keeptime}.}
#' \item{summary}{\code{list} containing a collection of
#' summary statistics of the posterior draws for \code{para}, \code{latent},
#' and \code{latent0}.}
#' \item{meanmodel}{\code{character} containing information about how \code{designmatrix}
#' was employed.}
#'
#' To display the output, use \code{print}, \code{summary} and \code{plot}. The
#' \code{print} method simply prints the posterior draws (which is very likely
#' a lot of output); the \code{summary} method displays the summary statistics
#' currently stored in the object; the \code{plot} method
#' \code{\link{plot.svdraws}} gives a graphical overview of the posterior
#' distribution by calling \code{\link{volplot}}, \code{\link{traceplot}} and
#' \code{\link{densplot}} and displaying the results on a single page.
#' @note If \code{y} contains zeros, you might want to consider de-meaning your
#' returns or use \code{designmatrix = "ar0"}.
#' @author Gregor Kastner \email{[email protected]@wu.ac.at}
#' @seealso \code{\link{svsim}}, \code{\link{svlsample}}, \code{\link{updatesummary}},
#' \code{\link{predict.svdraws}}, \code{\link{plot.svdraws}}.
#' @references Kastner, G. and Frühwirth-Schnatter, S. (2014).
#' Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC
#' estimation of stochastic volatility models. \emph{Computational Statistics &
#' Data Analysis}, \bold{76}, 408--423,
#' \url{http://dx.doi.org/10.1016/j.csda.2013.01.002}.
#' @keywords models ts
#' @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$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) #' #' #' \dontrun{ #' # 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) #' } #' @export svsample <- function(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, ...) { # Some error checking for y if (inherits(y, "svsim")) { y <- y[["y"]] warning("Extracted data vector from 'svsim'-object.") } if (!is.numeric(y)) stop("Argument 'y' (data vector) must be numeric.") if (length(y) < 2) stop("Argument 'y' (data vector) must contain at least two elements.") y_orig <- y y <- as.vector(y) if (any(is.na(designmatrix)) && any(y^2 == 0)) { myoffset <- sd(y)/10000 warning(paste("Argument 'y' (data vector) contains zeros. I am adding an offset constant of size ", myoffset, " to do the auxiliary mixture sampling. If you want to avoid this, you might consider de-meaning the returns before calling this function.", sep="")) } else myoffset <- 0 # Some error checking for draws if (!is.numeric(draws) || length(draws) != 1 || draws < 1) { stop("Argument 'draws' (number of MCMC iterations after burn-in) must be a single number >= 1.") } else { draws <- as.integer(draws) } # Some error checking for burnin if (!is.numeric(burnin) || length(burnin) != 1 || burnin < 0) { stop("Argument 'burnin' (burn-in period) must be a single number >= 0.") } else { burnin <- as.integer(burnin) } # Some error checking for designmatrix and save the mode of mean modelling meanmodel <- "matrix" arorder <- 0L if (any(is.na(designmatrix))) { designmatrix <- matrix(NA) meanmodel <- "none" } else { if (any(grep("ar[0-9]+$", as.character(designmatrix)[1]))) {
arorder <- as.integer(gsub("ar", "", as.character(designmatrix)))
if (length(y) <= arorder + 1) stop("Time series 'y' is too short for this AR process.")
designmatrix <- matrix(rep(1, length(y) - arorder), ncol = 1)
colnames(designmatrix) <- c("const")
meanmodel <- "constant"
if (arorder >= 1) {  # e.g. first row with "ar(3)": c(1, y_3, y_2, y_1)
for (i in 1:arorder) {
oldnames <- colnames(designmatrix)
designmatrix <- cbind(designmatrix, y[(arorder-i+1):(length(y)-i)])
colnames(designmatrix) <- c(oldnames, paste0("ar", i))
}
y <- y[-(1:arorder)]
meanmodel <- paste0("ar", arorder)
}
}
if (!is.numeric(designmatrix)) stop("Argument 'designmatrix' must be a numeric matrix or an AR-specification.")
if (!is.matrix(designmatrix)) {
designmatrix <- matrix(designmatrix, ncol = 1)
}
if (!nrow(designmatrix) == length(y)) stop("Number of columns of argument 'designmatrix' must be equal to length(y).")
}

# Some error checking for the prior parameters
if (!is.numeric(priormu) || length(priormu) != 2) {
stop("Argument 'priormu' (mean and sd for the Gaussian prior for mu) must be numeric and of length 2.")
}

if (!is.numeric(priorphi) | length(priorphi) != 2) {
stop("Argument 'priorphi' (shape1 and shape2 parameters for the Beta prior for (phi + 1) / 2) must be numeric and of length 2.")
}

if (!is.numeric(priorsigma) | length(priorsigma) != 1 | priorsigma <= 0) {
stop("Argument 'priorsigma' (scaling of the chi-squared(df = 1) prior for sigma^2) must be a single number > 0.")
}

if (any(is.na(priornu))) {
priornu <- NA
} else {
if (!is.numeric(priornu) || length(priornu) != 2 || priornu[1] >= priornu[2] || priornu[1] < 0) {
stop("If not NA, argument 'priornu' (lower and upper bounds for the uniform prior for the df) must be numeric and of length 2. Moreover, 0 <= priornu[1] < priornu[2].")
}
}

if (!is.numeric(priorbeta) || length(priorbeta) != 2) {
stop("Argument 'priorbeta' (means and sds for the independent Gaussian priors for beta) must be numeric and of length 2.")
}

if (!is.numeric(priorlatent0) || length(priorlatent0) != 1 || priorlatent0 < 0) {
if (priorlatent0 == "stationary") priorlatent0 <- -1L else
stop("Argument 'priorlatent0' must be 'stationary' or a single non-negative number.")
}

# Some error checking for thinpara
if (!is.numeric(thinpara) || length(thinpara) != 1 || thinpara < 1) {
stop("Argument 'thinpara' (thinning parameter for mu, phi, and sigma) must be a single number >= 1.")
} else {
thinpara <- as.integer(thinpara)
}

# Some error checking for thinlatent
if (!is.numeric(thinlatent) || length(thinlatent) != 1 || thinlatent < 1) {
stop("Argument 'thinlatent' (thinning parameter for the latent log-volatilities) must be a single number >= 1.")
} else {
thinlatent <- as.integer(thinlatent)
}

# Check whether 'thintime' was used
if (!is.null(thintime))
stop("Argument 'thintime' is deprecated. Please use 'keeptime' instead.")

# Some error checking for keeptime
if (length(keeptime) != 1L || !is.character(keeptime) || !(keeptime %in% c("all", "last"))) {
stop("Parameter 'keeptime' must be either 'all' or 'last'.")
} else {
if (keeptime == "all") thintime <- 1L else if (keeptime == "last") thintime <- length(y)
}

# Some error checking for expert
if (missing(expert)) {
para <- 3L ; parameterization <- 'GIS_C'
mhcontrol <- -1
gammaprior <- TRUE
truncnormal <- FALSE
mhsteps <- 2L
B011 <- 10^8
B022 <- 10^12
} else {
expertnames <- names(expert)
if (!is.list(expert) | is.null(expertnames) | any(expertnames == ""))
stop("Argument 'expert' must be a named list with nonempty names.")
if (length(unique(expertnames)) != length(expertnames))
stop("No duplicate elements allowed in argument 'expert'.")
allowednames <- c("parameterization", "mhcontrol", "gammaprior", "truncnormal", "mhsteps", "proposalvar4sigmaphi", "proposalvar4sigmatheta")
exist <- pmatch(expertnames, allowednames)
if (any(is.na(exist)))
stop(paste("Illegal element '", paste(expertnames[is.na(exist)], collapse="' and '"), "' in argument 'expert'.", sep=''))

expertenv <- list2env(expert)

if (exists("parameterization", expertenv)) {
parameterization <- expert[["parameterization"]]
if (!is.character(parameterization) || any(is.na(parameterization))) {
stop("Argument 'parameterization' must be either 'centered', 'noncentered', 'GIS_C', or 'GIS_NC'.")
}
switch(parameterization,
centered = para <- 1L,
noncentered = para <- 2L,
GIS_C = para <- 3L,
GIS_NC = para <- 4L,
stop("Unknown parameterization. Currently you can only use 'centered', 'noncentered', 'GIS_C', and 'GIS_NC'.")
)
} else {
para <- 3L ; parameterization <- 'GIS_C'
}

# Remark: mhcontrol < 0 means independence proposal,
#         mhcontrol > 0 controls stepsize of log-random-walk proposal
if (exists("mhcontrol", expertenv)) {
mhcontrol <- expert[["mhcontrol"]]
if (!is.numeric(mhcontrol) | length(mhcontrol) != 1)
stop("Argument 'mhcontrol' must be a single number.")
} else {
mhcontrol <- -1
}

# use a Gamma prior for sigma^2 in C?
if (exists("gammaprior", expertenv)) {
gammaprior <- expert[["gammaprior"]]
if (!is.logical(gammaprior)) stop("Argument 'gammaprior' must be TRUE or FALSE.")
} else {
gammaprior <- TRUE
}

# use a truncated normal as proposal? (or normal with rejection step)
if (exists("truncnormal", expertenv)) {
truncnormal <- expert[["truncnormal"]]
if (!is.logical(truncnormal)) stop("Argument 'truncnormal' must be TRUE or FALSE.")
} else {
truncnormal <- FALSE
}

if (exists("mhsteps", expertenv)) {
mhsteps <- as.integer(expert[["mhsteps"]])
if (mhsteps != 2L & mhsteps != 1L & mhsteps != 3L) stop("mhsteps must be 1, 2, or 3")
if (mhsteps != 2L & mhcontrol >= 0)
stop("Log normal random walk proposal currently only implemented for mhsteps==2.")
if (mhsteps != 2L & !isTRUE(gammaprior))
stop("Inverse Gamma prior currently only implemented for mhsteps==2.")
} else {
mhsteps <- 2L
}

# prior for ridge _proposal_ (variance of sigma*phi)
if (exists("proposalvar4sigmaphi", expertenv)) {
B011 <- expert[["proposalvar4sigmaphi"]]
if (!is.numeric(B011) | length(B011) != 1 | B011 <= 0)
stop("Argument 'proposalvar4sigmaphi' must be a positive number.")
} else {
B011 <- 10^8
}

# prior for ridge _proposal_ (variance of sigma*theta)
if (exists("proposalvar4sigmatheta", expertenv)) {
B022 <- expert[["proposalvar4sigmatheta"]]
if (!is.numeric(B022) | length(B022) != 1 | B022 <= 0)
stop("Argument 'proposalvar4sigmatheta' must be a positive number.")
} else {
B022 <- 10^12
}
}

# Some input checking for startpara
if (missing(startpara)) {
if (any(is.na(priornu))) {
startpara <- list(mu = priormu[1],
phi = 2 * (priorphi[1] / sum(priorphi)) - 1,
sigma = priorsigma)
} else {
startpara <- list(mu = priormu[1],
phi = 2 * (priorphi[1] / sum(priorphi)) - 1,
sigma = priorsigma,
nu = mean(priornu))
}
} else {
if (!is.list(startpara))
stop("Argument 'startpara' must be a list. Its elements must be named 'mu', 'phi', 'sigma'. Moreover, if !is.na(priornu), an element named 'nu' must exist.")

if (!is.numeric(startpara[["mu"]]))
stop('Argument \'startpara[["mu"]]\' must exist and be numeric.')

if (!is.numeric(startpara[["phi"]]))
stop('Argument \'startpara[["phi"]]\' must exist and be numeric.')

if (abs(startpara[["phi"]]) >= 1)
stop('Argument \'startpara[["phi"]]\' must be between -1 and 1.')

if (!is.numeric(startpara[["sigma"]]))
stop('Argument \'startpara[["sigma"]]\' must exist and be numeric.')

if (startpara[["sigma"]] <= 0)
stop('Argument \'startpara[["sigma"]]\' must be positive.')

if (!any(is.na(priornu)) && !is.numeric(startpara[["nu"]]))
stop('Argument \'startpara[["nu"]]\' must exist and be numeric.')

if (!any(is.na(priornu)) && (startpara[["nu"]] > priornu[2] || startpara[["nu"]] < priornu[1]))
stop('Argument \'startpara[["nu"]]\' must be within range(priornu).')
}

# Some input checking for startlatent
if (missing(startlatent)) {
startlatent <- rep(-10, length(y))
} else {
if (!is.numeric(startlatent) | length(startlatent) != length(y))
stop("Argument 'startlatent' must be numeric and of the same length as the data 'y'.")
}

if (length(keeptau) != 1 || !is.logical(keeptau)) {
stop("Argument 'keeptau' must be TRUE or FALSE.")
}

if (any(is.na(priornu)) && keeptau) {
warning("Setting argument 'keeptau' to FALSE, as 'priornu' is NA.")
keeptau <- FALSE
}

if (!quiet) {
cat(paste("\nCalling ", parameterization, " MCMC sampler with ", draws+burnin, " iter. Series length is ", length(y), ".\n",sep=""), file=stderr())
flush.console()
}

if (.Platform$OS.type != "unix") myquiet <- TRUE else myquiet <- quiet # Hack to prevent console flushing problems with Windows runtime <- system.time(res <- svsample_cpp(y, draws, burnin, designmatrix, priormu[1], priormu[2]^2, priorphi[1], priorphi[2], priorsigma, thinlatent, thintime, startpara, startlatent, keeptau, myquiet, para, mhsteps, B011, B022, mhcontrol, gammaprior, truncnormal, myoffset, FALSE, priornu, priorbeta, priorlatent0)) if (any(is.na(res))) stop("Sampler returned NA. This is most likely due to bad input checks and shouldn't happen. Please report to package maintainer.") if (!quiet) { cat("Timing (elapsed): ", file=stderr()) cat(runtime["elapsed"], file=stderr()) cat(" seconds.\n", file=stderr()) cat(round((draws+burnin)/runtime[3]), "iterations per second.\n\n", file=stderr()) cat("Converting results to coda objects... ", file=stderr()) } # store results: # remark: +1, because C-sampler also returns the first value res$y <- y_orig
res$para <- mcmc(res$para[seq(burnin+thinpara+1, burnin+draws+1, thinpara),,drop=FALSE], burnin+thinpara, burnin+draws, thinpara)
res$latent <- mcmc(t(res$latent), burnin+thinlatent, burnin+draws, thinlatent)
if (keeptime == "all") attr(res$latent, "dimnames") <- list(NULL, paste0('h_', arorder + seq_along(y))) else if (keeptime == "last") attr(res$latent, "dimnames") <- list(NULL, paste0('h_', arorder + length(y)))
res$latent0 <- mcmc(res$latent0, burnin+thinlatent, burnin+draws, thinlatent)
if (!any(is.na(designmatrix))) {
res$beta <- mcmc(res$beta[seq(burnin+thinpara+1, burnin+draws+1, thinpara),,drop=FALSE], burnin+thinpara, burnin+draws, thinpara)
attr(res$beta, "dimnames") <- list(NULL, paste("b", 0:(ncol(designmatrix)-1), sep = "_")) } else res$beta <- NULL
res$meanmodel <- meanmodel if (ncol(res$para) == 3) {
attr(res$para, "dimnames") <- list(NULL, c("mu", "phi", "sigma")) res$priors <- list(mu = priormu, phi = priorphi, sigma = priorsigma, gammaprior = gammaprior)
} else {
attr(res$para, "dimnames") <- list(NULL, c("mu", "phi", "sigma", "nu")) res$priors <- list(mu = priormu, phi = priorphi, sigma = priorsigma, nu = priornu, gammaprior = gammaprior)
#res$tau <- mcmc(t(res$tau), burnin+thinlatent, burnin+draws, thinlatent)
}

if (!any(is.na(designmatrix))) {
res$priors <- c(res$priors, "beta" = list(priorbeta), "designmatrix" = list(designmatrix))
}

if (keeptau) {
res$tau <- mcmc(t(res$tau), burnin+thinlatent, burnin+draws, thinlatent)
if (keeptime == "all") attr(res$tau, "dimnames") <- list(NULL, paste0('tau_', arorder + seq_along(y))) else if (keeptime == "last") attr(res$tau, "dimnames") <- list(NULL, paste0('tau_', arorder + length(y)))
}

res$runtime <- runtime res$thinning <- list(para = thinpara, latent = thinlatent, time = keeptime)
class(res) <- "svdraws"

if (!quiet) {
cat("Done!\n", file=stderr())
cat("Summarizing posterior draws... ", file=stderr())
}
res <- updatesummary(res, ...)

if (!quiet) cat("Done!\n\n", file=stderr())
res
}

# This function does not check input nor converts the result to coda objects!

#' @rdname svsample
#' @export
svtsample <- function(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, ...) {
svsample(y, draws = draws, burnin = burnin, designmatrix = designmatrix,
priormu = priormu, priorphi = priorphi, priorsigma = priorsigma,
priornu = priornu, priorbeta = priorbeta, priorlatent0 = priorlatent0,
thinpara = thinpara, thinlatent = thinlatent, keeptime = keeptime, thintime = thintime,
keeptau = keeptau, quiet = quiet, startpara = startpara, startlatent = startlatent,
expert = expert, ...)
}

#' Minimal overhead version of \code{\link{svsample}}.
#'
#' \code{svsample2} is a minimal overhead version of \code{\link{svsample}}
#' with slightly different default arguments and a simplified return value
#' structure. It is intended to be used mainly for one-step updates where speed
#' is an issue, e.g., as a plug-in into other MCMC samplers. Note that
#' absolutely no input checking is performed, thus this function is to be used
#' with proper care!
#'
#' As opposed to the ordinary \code{\link{svsample}}, the default values differ
#' for \code{draws}, \code{burnin}, and \code{quiet}. Note that currently
#' neither \code{expert} nor \code{\dots{}} arguments are provided.
#'
#' @param y numeric vector containing the data (usually log-returns), which
#' must not contain zeroes.
#' @param 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 defaults value is 1.
#' @param 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 0.
#' @param priormu numeric vector of length 2, indicating mean and standard
#' deviation for the Gaussian prior distribution of the parameter \code{mu},
#' the level of the log-volatility. The default value is \code{c(0, 100)},
#' which constitutes a practically uninformative prior for common exchange rate
#' datasets, stock returns and the like.
#' @param priorphi numeric vector of length 2, indicating the shape parameters
#' for the Beta prior distribution of the transformed parameter
#' \code{(phi + 1) / 2}, where \code{phi} denotes the persistence of the
#' log-volatility. The default value is \code{c(5, 1.5)}, which constitutes a
#' prior that puts some belief in a persistent log-volatility but also
#' encompasses the region where \code{phi} is around 0.
#' @param priorsigma single positive real number, which stands for the scaling
#' of the transformed parameter \code{sigma^2}, where \code{sigma} denotes the
#' volatility of log-volatility. More precisely, \code{sigma^2 ~ priorsigma *
#' chisq(df = 1)}. The default value is \code{1}, which constitutes a
#' reasonably vague prior for many common exchange rate datasets, stock returns
#' and the like.
#' @param priornu numeric vector of length 2 (or \code{NA}), indicating the
#' lower and upper bounds for the uniform prior distribution of the parameter
#' \code{nu}, the degrees-of-freedom parameter of the conditional innovations
#' t-distribution. The default value is \code{NA}, fixing the
#' degrees-of-freedom to infinity. This corresponds to conditional standard
#' normal innovations, the pre-1.1.0 behavior of \pkg{stochvol}.
#' @param priorlatent0 either a single non-negative number or the string
#' \code{'stationary'} (the default, also the behavior before version 1.3.0).
#' When \code{priorlatent0} is equal to \code{'stationary'}, the stationary
#' distribution of the latent AR(1)-process is used as the prior for the
#' initial log-volatility \code{h_0}. When \code{priorlatent0} is equal to a
#' number \eqn{B}, we have \eqn{h_0 \sim N(\mu, B\sigma^2)} a priori.
#' @param thinpara single number greater or equal to 1, coercible to integer.
#' Every \code{thinpara}th parameter draw is kept and returned. The default
#' value is 1, corresponding to no thinning of the parameter draws -- every
#' draw is stored.
#' @param thinlatent single number greater or equal to 1, coercible to integer.
#' Every \code{thinlatent}th 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.
#' @param thintime \emph{Deprecated.} Use 'keeptime' instead.
#' @param keeptime Either 'all' (the default) or 'last'. Indicates which latent
#  volatility draws should be stored.
#' @param 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'.
#' @param quiet logical value indicating whether the progress bar and other
#' informative output during sampling should be omitted. The default value is
#' \code{TRUE}, implying non-verbose output.
#' @param startpara \emph{compulsory} named list, containing the starting
#' values for the parameter draws. \code{startpara} must contain three elements
#' named \code{mu}, \code{phi}, and \code{sigma}, where \code{mu} is an
#' arbitrary numerical value, \code{phi} is a real number between \code{-1} and
#' \code{1}, and \code{sigma} is a positive real number. Moreover, if
#' \code{priornu} is not \code{NA}, \code{startpara} must also contain an
#' element named \code{nu} (the degrees of freedom parameter for the
#' t-innovations).
#' @param startlatent \emph{compulsory} vector of length \code{length(x$y)}, #' containing the starting values for the latent log-volatility draws. #' @return A list with three components: #' \item{para}{\code{3} times #' \code{draws} matrix containing the parameter draws. If \code{priornu} is not #' \code{NA}, this is a \code{4} times \code{draws} matrix.} #' \item{latent}{\code{length(y)} times \code{draws} matrix containing draws of #' the latent variables \code{h_1, \dots{}, h_n}.} #' \item{latent0}{Vector of #' length \code{draws} containing the draw(s) of the initial latent variable #' \code{h_0}.} #' @note Please refer to the package vignette for an example. #' @section Warning: Expert use only! For most applications, the use of #' \code{\link{svsample}} is recommended. #' @author Gregor Kastner \email{[email protected]@wu.ac.at} #' @seealso \code{\link{svsample}} #' @keywords models ts #' @examples #' data(exrates) #' aud.price <- subset(exrates, #' as.Date("2010-01-01") <= date & date < as.Date("2011-01-01"), #' "AUD")[,1] #' draws <- svsample2(logret(aud.price), #' draws = 10, burnin = 0, #' startpara = list(phi = 0.95, mu = -10, sigma = 0.2, rho = -0.1), #' startlatent = rep_len(-10, length(aud.price) - 1)) #' @export svsample2 <- function(y, draws = 1, burnin = 0, priormu = c(0, 100), priorphi = c(5, 1.5), priorsigma = 1, priornu = NA, priorlatent0 = "stationary", thinpara = 1, thinlatent = 1, thintime = NULL, keeptime = "all", keeptau = FALSE, quiet = TRUE, startpara, startlatent) { if (priorlatent0 == "stationary") priorlatent0 <- -1L # Check whether 'thintime' was used if (!is.null(thintime)) stop("Argument 'thintime' is deprecated. Please use 'keeptime' instead.") if (keeptime == "all") thintime <- 1L else if (keeptime == "last") thintime <- length(y) else stop("Parameter 'keeptime' must be either 'all' or 'last'.") res <- svsample_cpp(y, draws, burnin, matrix(NA), priormu[1], priormu[2]^2, priorphi[1], priorphi[2], priorsigma, thinlatent, thintime, startpara, startlatent, keeptau, quiet, 3L, 2L, 10^8, 10^12, -1, TRUE, FALSE, 0, FALSE, priornu, c(NA, NA), priorlatent0) res$para <- t(res$para[seq(burnin+thinpara+1, burnin+draws+1, thinpara),,drop=FALSE]) if (nrow(res$para) == 3) {
rownames(res$para) <- names(res$para) <- c("mu", "phi", "sigma")
} else {
rownames(res$para) <- names(res$para) <- c("mu", "phi", "sigma", "nu")
}
res.meanmodel <- "none"

res
}

#' Markov Chain Monte Carlo (MCMC) Sampling for the Stochastic Volatility
#' Model with Leverage (SVL)
#'
#' \code{svlsample} simulates from the joint posterior distribution of the SVL
#' parameters \code{mu}, \code{phi}, \code{sigma}, and \code{rho},
#' along with the latent log-volatilities \code{h_1,...,h_n} and returns the
#' MCMC draws. If a design matrix is provided, simple Bayesian regression can
#' also be conducted.
#'
#' @param y numeric vector containing the data (usually log-returns), which
#' must not contain zeros. Alternatively, \code{y} can be an \code{svsim}
#' object. In this case, the returns will be extracted and a warning is thrown.
#' @param 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.
#' @param 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.
#' @param designmatrix regression design matrix for modeling the mean. Must
#' have \code{length(y)} rows. Alternatively, \code{designmatrix} may be a
#' string of the form \code{"arX"}, where \code{X} is a nonnegative integer. To
#' fit a constant mean model, use \code{designmatrix = "ar0"} (which is
#' equivalent to \code{designmatrix = matrix(1, nrow = length(y))}). To fit an
#' AR(1) model, use \code{designmatrix = "ar1"}, and so on. If some elements of
#' \code{designmatrix} are \code{NA}, the mean is fixed to zero (pre-1.2.0
#' behavior of \pkg{stochvol}).
#' @param priormu numeric vector of length 2, indicating mean and standard
#' deviation for the Gaussian prior distribution of the parameter \code{mu},
#' the level of the log-volatility. The default value is \code{c(0, 100)},
#' which constitutes a practically uninformative prior for common exchange rate
#' datasets, stock returns and the like.
#' @param priorphi numeric vector of length 2, indicating the shape parameters
#' for the Beta prior distribution of the transformed parameter
#' \code{(phi + 1) / 2}, where \code{phi} denotes the persistence of the
#' log-volatility. The default value is \code{c(5, 1.5)}, which constitutes a
#' prior that puts some belief in a persistent log-volatility but also
#' encompasses the region where \code{phi} is around 0.
#' @param priorsigma single positive real number, which stands for the scaling
#' of the transformed parameter \code{sigma^2}, where \code{sigma} denotes the
#' volatility of log-volatility. More precisely, \code{sigma^2 ~ priorsigma *
#' chisq(df = 1)}. The default value is \code{1}, which constitutes a
#' reasonably vague prior for many common exchange rate datasets, stock returns
#' and the like.
#' @param priorrho numeric vector of length 2, indicating the shape parameters
#' for the Beta prior distribution of the transformed parameter
#' \code{(rho + 1) / 2}, where \code{rho} denotes the conditional correlation
#' between observation and the increment of the
#' log-volatility. The default value is \code{c(4, 4)}, which constitutes a
#' slightly informative prior around 0 (the no leverage case) to boost convergence.
#' @param priorbeta numeric vector of length 2, indicating the mean and
#' standard deviation of the Gaussian prior for the regression parameters. The
#' default value is \code{c(0, 10000)}, which constitutes a very vague prior
#' for many common datasets. Not used if \code{designmatrix} is \code{NA}.
#' @param thinpara single number greater or equal to 1, coercible to integer.
#' Every \code{thinpara}th 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.
#' @param thinlatent single number greater or equal to 1, coercible to integer.
#' Every \code{thinlatent}th 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.
#' @param thintime \emph{Deprecated.} Use 'keeptime' instead.
#' @param keeptime Either 'all' (the default) or 'last'. Indicates which latent
#  volatility draws should be stored.
#' @param quiet logical value indicating whether the progress bar and other
#' informative output during sampling should be omitted. The default value is
#' \code{FALSE}, implying verbose output.
#' @param startpara \emph{optional} named list, containing the starting values
#' for the parameter draws. If supplied, \code{startpara} must contain four
#' elements named \code{mu}, \code{phi}, \code{sigma}, and \code{rho}, where \code{mu} is
#' an arbitrary numerical value, \code{phi} is a real number between \code{-1}
#' and \code{1}, \code{sigma} is a positive real number, and \code{rho} is
#' a real number between \code{-1} and \code{1}. The default value is equal
#' to the prior mean.
#' @param startlatent \emph{optional} vector of length \code{length(y)},
#' containing the starting values for the latent log-volatility draws. The
#' default value is \code{rep(-10, length(y))}.
#' @param expert \emph{optional} named list of expert parameters. For most
#' applications, the default values probably work best. If
#' \code{expert} is provided, it may contain the following named elements:
#'
#' \code{parameterization}: Character string containing values \code{"centered"},
#' and \code{"noncentered"}. Alternatively, it can be a single element character
#' vector of the form \code{"asisX"}, where \code{X} is an integer, which is
#' equivalent to \code{rep(c("centered", "noncentered"), X)}.
#' Defaults to \code{"asis5"}.
#'
#' \code{gammaprior}: Single logical value indicating whether a Gamma prior for
#' \code{sigma^2} should be used. If set to \code{FALSE}, a moment-matched Inverse Gamma
#' prior is employed. Defaults to \code{TRUE}.
#'
#' \code{init.with.svsample}: Single integer indicating the length of a pre-burnin'' run using
#' the computationally much more efficient \code{\link{svsample}}. This run helps
#' in finding good initial values for the latent states, giving \code{svlsample}
#' a considerable initial boost for convergence. Defaults to \code{1000L}.
#'
#' \code{mhcontrol}: Either a single numeric value specifying the diagonal elements of
#' a diagonal covariance matrix, or a list with two elements, both single numeric values
#' (explained later), or a 4x4 covariance matrix. Argument \code{mhcontrol} controls the
#' proposal density of a Metropolis-Hastings (MH) update step when jointly sampling \code{mu},
#' \code{phi}, \code{sigma}, and \code{rho}. It specifies the covariance matrix of a
#' log-random-walk proposal. In case \code{mhcontrol} is a list of length two, its elements
#' have to be \code{scale} and \code{rho.var}. In this case, the covariance matrix is calculated
#' from the pre-burnin step with \code{\link{svsample}}, which gives an approximate
#' posterior structure of the second moment for \code{mu}, \code{phi}, and \code{sigma}.
#' This covariance matrix is then extended with \code{mhcontrol$rho.var}, specifying the #' variance for \code{rho}. The off-diagonal elements belonging to \code{rho} are set #' to 0. Finally, the whole covariance matrix is scaled by \code{mhcontrol$scale}. For
#' this case to work, \code{init.with.svsample} has to be positive.
#' Defaults to \code{list(scale=0.35, rho.var=0.02)}.
#'
#' \code{correct.latent.draws}: Single logical value indicating whether to correct
#' the draws obtained from the auxiliary model of Omori, et al. (2007). Defaults
#' to \code{TRUE}.
#' @param \dots Any extra arguments will be forwarded to
#' \code{\link{updatesummary}}, controlling the type of statistics calculated
#' for the posterior draws.
#' @return The value returned is a list object of class \code{svldraws} holding
#' \item{para}{\code{mcmc} object containing the \emph{parameter} draws from
#' the posterior distribution.}
#' \item{latent}{\code{mcmc} object containing the
#' \emph{latent instantaneous log-volatility} draws from the posterior
#' distribution.}
#' \item{beta}{\code{mcmc} object containing the \emph{regression coefficient}
#' draws from the posterior distribution \emph{(optional)}.}
#' \item{y}{the argument \code{y}.}
#' \item{runtime}{\code{proc_time} object containing the
#' run time of the sampler.}
#' \item{priors}{\code{list} containing the parameter
#' values of the prior distribution, i.e. the arguments \code{priormu},
#' \code{priorphi}, \code{priorsigma}, and \code{priorrho}, and potentially
#' \code{priorbeta}.}
#' \item{thinning}{\code{list} containing the thinning
#' parameters, i.e. the arguments \code{thinpara}, \code{thinlatent} and
#' \code{keeptime}.}
#' \item{summary}{\code{list} containing a collection of
#' summary statistics of the posterior draws for \code{para}, and \code{latent}.}
#' \item{meanmodel}{\code{character} containing information about how \code{designmatrix}
#' was used.}
#'
#' To display the output, use \code{print}, \code{summary} and \code{plot}. The
#' \code{print} method simply prints the posterior draws (which is very likely
#' a lot of output); the \code{summary} method displays the summary statistics
#' currently stored in the object; the \code{plot} method
#' \code{\link{plot.svdraws}} gives a graphical overview of the posterior
#' distribution by calling \code{\link{volplot}}, \code{\link{traceplot}} and
#' \code{\link{densplot}} and displaying the results on a single page.
#' @note If \code{y} contains zeros, you might want to consider de-meaning your
#' returns or use \code{designmatrix = "ar0"}. We use the Metropolis-Hastings
#' algorithm for sampling the latent vector \code{h}, where the proposal is a
#' draw from an auxiliary mixture approximation model [Omori, et al. (2007)].
#' We draw the parameters \code{mu}, \code{phi}, \code{sigma}, and \code{rho}
#' jointly by employing a Metropolis random walk step. By default, we boost the
#' random walk through the repeated application of the ancillarity-sufficiency
#' interweaving strategy (ASIS) [Yu, Meng (2011)]. A message in the beginning
#' of sampling indicates the interweaving strategy used, which can be modified
#' through parameter \code{expert}.
#' @author Darjus Hosszejni \email{[email protected]@wu.ac.at}
#' @references
#' Yu, Y. and Meng, X.-L. (2011).
#' To Center or not to Center: That is not the Question---An Ancillarity-Sufficiency
#' Interweaving Strategy (ASIS) for Boosting MCMC Efficiency. \emph{Journal of
#' Computational and Graphical Statistics}, \bold{20}(3), 531--570,
#' \url{http://dx.doi.org/10.1198/jcgs.2011.203main}
#'
#' Omori, Y. and Chib, S. and Shephard, N. and Nakajima, J. (2007).
#' Stochastic Volatility with Leverage: Fast and Efficient Likelihood Inference.
#' \emph{Journal of Econometrics}, \bold{140}(2), 425--449,
#' \url{http://dx.doi.org/10.1016/j.jeconom.2006.07.008}
#' @seealso \code{\link{svsim}}, \code{\link{svsample}}, \code{\link{updatesummary}},
#' \code{\link{predict.svdraws}}, \code{\link{plot.svdraws}}.
#' @keywords models ts
#' @examples
#' \dontrun{
#' # Example 1
#' ## Simulate a short SVL process
#' sim <- svsim(200, mu = -10, phi = 0.95, sigma = 0.2, rho = -0.4)
#'
#' ## Obtain 5000 draws from the sampler (that's not a lot)
#' draws <- svlsample(sim$y) #' #' ## Check out the results #' summary(draws) #' plot(draws, simobj = sim) #' #' #' # Example 2 #' ## AR(1) structure for the mean #' data(exrates) #' len <- 1200 #' ahead <- 100 #' y <- head(exrates$USD, len)
#'
#' ## Fit AR(1)-SVL model to EUR-USD exchange rates
#' res <- svlsample(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 <- 1200
#' 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 <- svlsample(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) #' } #' @export svlsample <- function (y, draws = 10000, burnin = 1000, designmatrix = NA, priormu = c(0, 100), priorphi = c(5, 1.5), priorsigma = 1, priorrho = c(4, 4), priorbeta = c(0, 10000), thinpara = 1, thinlatent = 1, thintime = NULL, keeptime = "all", quiet = FALSE, startpara, startlatent, expert, ...) { # Some error checking for y if (inherits(y, "svsim")) { y <- y[["y"]] warning("Extracted data vector from 'svsim'-object.") } if (!is.numeric(y)) stop("Argument 'y' (data vector) must be numeric.") if (length(y) < 2) stop("Argument 'y' (data vector) must contain at least two elements.") y_orig <- y y <- as.vector(y) myoffset <- if (any(is.na(designmatrix)) && any(y^2 == 0)) 1.0e-8 else 0 if (myoffset > 0) { warning(paste("Argument 'y' (data vector) contains zeros. I am adding an offset constant of size ", myoffset, " to do the auxiliary mixture sampling. If you want to avoid this, you might consider de-meaning the returns before calling this function.", sep="")) } # Some error checking for draws if (!is.numeric(draws) || length(draws) != 1 || draws < 1) { stop("Argument 'draws' (number of MCMC iterations after burn-in) must be a single number >= 1.") } else { draws <- as.integer(draws) } # Some error checking for burnin if (!is.numeric(burnin) || length(burnin) != 1 || burnin < 0) { stop("Argument 'burnin' (burn-in period) must be a single number >= 0.") } else { burnin <- as.integer(burnin) } # Some error checking for designmatrix meanmodel <- "matrix" arorder <- 0L if (any(is.na(designmatrix))) { designmatrix <- matrix(NA) meanmodel <- "none" } else { if (any(grep("ar[0-9]+$", as.character(designmatrix)[1]))) {
arorder <- as.integer(gsub("ar", "", as.character(designmatrix)))
if (length(y) <= (arorder + 1L)) stop("Time series 'y' is to short for this AR process.")
designmatrix <- matrix(rep(1, length(y) - arorder), ncol = 1)
colnames(designmatrix) <- c("const")
meanmodel <- "constant"
if (arorder >= 1) {
for (i in seq_len(arorder)) {
oldnames <- colnames(designmatrix)
designmatrix <- cbind(designmatrix, y[(arorder-i+1):(length(y)-i)])
colnames(designmatrix) <- c(oldnames, paste0("ar", i))
}
y <- y[-(1:arorder)]
meanmodel <- paste0("ar", arorder)
}
}
if (!is.numeric(designmatrix)) stop("Argument 'designmatrix' must be a numeric matrix or an AR-specification.")
if (!is.matrix(designmatrix)) {
designmatrix <- matrix(designmatrix, ncol = 1)
}
if (!nrow(designmatrix) == length(y)) stop("Number of columns of argument 'designmatrix' must be equal to length(y).")
}

# Some error checking for the prior parameters
if (!is.numeric(priormu) || length(priormu) != 2) {
stop("Argument 'priormu' (mean and variance for the Gaussian prior for mu) must be numeric and of length 2.")
}

if (!is.numeric(priorphi) || length(priorphi) != 2) {
stop("Argument 'priorphi' (shape1 and shape2 parameters for the Beta prior for (phi + 1) / 2) must be numeric and of length 2.")
}

if (!is.numeric(priorsigma) || length(priorsigma) != 1 || priorsigma <= 0) {
stop("Argument 'priorsigma' (scaling of the chi-squared(df = 1) prior for sigma^2) must be a single number > 0.")
}

if (!is.numeric(priorrho) || length(priorrho) != 2) {
stop("Argument 'priorrho' (shape1 and shape2 parameters for the Beta prior for (rho + 1) / 2) must be numeric and of length 2.")
}

if (!is.numeric(priorbeta) || length(priorbeta) != 2) {
stop("Argument 'priorbeta' (means and sds for the independent Gaussian priors for beta) must be numeric and of length 2.")
}

# Some error checking for thinpara
if (!is.numeric(thinpara) || length(thinpara) != 1 || thinpara < 1) {
stop("Argument 'thinpara' (thinning parameter for mu, phi, sigma, and rho) must be a single number >= 1.")
} else {
thinpara <- as.integer(thinpara)
}

# Some error checking for thinlatent
if (!is.numeric(thinlatent) || length(thinlatent) != 1 || thinlatent < 1) {
stop("Argument 'thinlatent' (thinning parameter for the latent log-volatilities) must be a single number >= 1.")
} else {
thinlatent <- as.integer(thinlatent)
}

# Check whether 'thintime' was used
if (!is.null(thintime))
stop("Argument 'thintime' is deprecated. Please use 'keeptime' instead.")

# Some error checking for keeptime
if (length(keeptime) != 1L || !is.character(keeptime) || !(keeptime %in% c("all", "last"))) {
stop("Parameter 'keeptime' must be either 'all' or 'last'.")
} else {
if (keeptime == "all") thintime <- 1L else if (keeptime == "last") thintime <- length(y)
}

# Some input checking for startpara
startparadefault <- list(mu = priormu[1],
phi = 2 * (priorphi[1] / sum(priorphi)) - 1,
sigma = priorsigma,
rho = 2 * (priorrho[1] / sum(priorrho)) - 1)
if (missing(startpara)) {
startpara <- startparadefault
} else {
if (!is.list(startpara))
stop("Argument 'startpara' must be a list. Its elements must be named 'mu', 'phi', 'sigma', and 'rho'.")

if (!is.numeric(startpara[["mu"]]))
stop('Argument \'startpara[["mu"]]\' must exist and be numeric.')

if (!is.numeric(startpara[["phi"]]))
stop('Argument \'startpara[["phi"]]\' must exist and be numeric.')

if (abs(startpara[["phi"]]) >= 1)
stop('Argument \'startpara[["phi"]]\' must be between -1 and 1.')

if (!is.numeric(startpara[["sigma"]]))
stop('Argument \'startpara[["sigma"]]\' must exist and be numeric.')

if (startpara[["sigma"]] <= 0)
stop('Argument \'startpara[["sigma"]]\' must be positive.')

if (!is.numeric(startpara[["rho"]]))
stop('Argument \'startpara[["rho"]]\' must exist and be numeric.')

if (abs(startpara[["rho"]]) >= 1)
stop('Argument \'startpara[["rho"]]\' must be between -1 and 1.')
}

# Some input checking for startlatent
if (missing(startlatent)) {
startlatent <- rep(-10, length(y))
} else {
if (!is.numeric(startlatent) | length(startlatent) != length(y))
stop("Argument 'startlatent' must be numeric and of the same length as the data 'y'.")
}

# Some error checking for expert
strategies <- c("centered", "noncentered")
expertdefault <- list(parameterization = rep(strategies, 5),  # default: ASISx5
mhcontrol = list(scale=0.35, rho.var=0.02),
gammaprior = TRUE,
init.with.svsample = 1000L,
correct.latent.draws = TRUE)
if (missing(expert)) {
parameterization <- expertdefault$parameterization mhcontrol <- expertdefault$mhcontrol
gammaprior <- expertdefault$gammaprior init.with.svsample <- expertdefault$init.with.svsample
correct.latent.draws <- expertdefault$correct.latent.draws } else { expertnames <- names(expert) if (!is.list(expert) || is.null(expertnames) || any(expertnames == "")) stop("Argument 'expert' must be a named list with nonempty names.") if (length(unique(expertnames)) != length(expertnames)) stop("No duplicate elements allowed in argument 'expert'.") allowednames <- c("parameterization", "mhcontrol", "gammaprior", "init.with.svsample", "correct.latent.draws") exist <- pmatch(expertnames, allowednames) if (any(is.na(exist))) stop(paste("Illegal element '", paste(expertnames[is.na(exist)], collapse="' and '"), "' in argument 'expert'.", sep='')) expertenv <- list2env(expert) if (exists("parameterization", expertenv)) { if (!is.character(expert[["parameterization"]])) stop("Argument 'parameterization' must be either a vector of 'centered', 'noncentered' values or a character string of form 'asis#' with # a positive integer.") nmatches <- grep("^asis[1-9][0-9]*$", expert[["parameterization"]])
if (length(nmatches) == 0) {
parameterization <- match.arg(expert[["parameterization"]], strategies, several.ok = TRUE)
} else if (length(nmatches) == 1) {
parameterization <- rep(strategies, nmatches)
} else {
parameterization <- NA
}
if (!all(parameterization %in% strategies)) {
stop("Argument 'parameterization' must be either a vector of 'centered', 'noncentered' values or a character string of form 'asis#' with # a positive integer.")
}
} else {
parameterization <- expertdefault$parameterization } if (exists("mhcontrol", expertenv)) { mhcontrol <- expert[["mhcontrol"]] if (!((is.numeric(mhcontrol) && ((length(mhcontrol) == 1 && mhcontrol > 0) || (is.matrix(mhcontrol) && NROW(mhcontrol) == 4 && NCOL(mhcontrol) == 4))) || (is.list(mhcontrol) && length(mhcontrol) == 2 && all(sort(c("row.var", "scale")) == sort(names(mhcontrol))) && is.numeric(mhcontrol$row.var) && mhcontrol$row.var > 0 && is.numeric(mhcontrol$scale) && mhcontrol$scale > 0))) stop("Argument 'mhcontrol' must be a single positive number, a list of two of positive numbers with names 'scale' and 'row.var', or a 4x4 covariance matrix.") } else { mhcontrol <- expertdefault$mhcontrol
}

if (exists("gammaprior", expertenv)) {
gammaprior <- expert[["gammaprior"]]
if (!is.logical(gammaprior)) stop("Argument 'gammaprior' must be TRUE or FALSE.")
} else {
gammaprior <- expertdefault$gammaprior } if (exists("init.with.svsample", expertenv)) { init.with.svsample <- expert[["init.with.svsample"]] if (!is.numeric(init.with.svsample) || length(init.with.svsample) != 1 || init.with.svsample < 0) stop("Argument 'init.with.svsample' must be a single non-negative number.") if (length(mhcontrol) == 2 && init.with.svsample == 0) stop("In case argument 'init.with.svsample' is set to 0, 'mhcontrol' cannot be of length 2.") init.with.svsample <- as.integer(init.with.svsample) } else { init.with.svsample <- expertdefault$init.with.svsample
}

if (exists("correct.latent.draws", expertenv)) {
correct.latent.draws <- expert[["correct.latent.draws"]]
if (!is.logical(correct.latent.draws)) stop("Argument 'correct.latent.draws' must be TRUE or FALSE.")
} else {
correct.latent.draws <- expertdefault$correct.latent.draws } } if (init.with.svsample > 0L) { if (!quiet) { cat(paste0("\nInitial values: calling 'svsample' with ", init.with.svsample, " iter. Series length is ", length(y), ".\n"), file=stderr()) flush.console() } init.runtime <- system.time({ init.res <- svsample(y, # not svsample2 because of 'designmatrix' draws = as.integer(init.with.svsample/2), burnin = as.integer(init.with.svsample/2), designmatrix = designmatrix, priormu = priormu, priorphi = priorphi, priorsigma = priorsigma, priornu = NA, priorbeta = priorbeta, priorlatent0 = "stationary", thinpara = 1L, thinlatent = 1L, keeptime = "all", keeptau = FALSE, quiet = TRUE, startpara = startpara[c("phi", "mu", "sigma")], startlatent = startlatent, expert = list(gammaprior = gammaprior, parameterization = "GIS_C"), quantiles = .5, esspara = FALSE, # fast updatesummary esslatent = FALSE) }) if (!quiet) { cat(paste0("Initial values: time taken by 'svsample': ", round(init.runtime["elapsed"], 3), " seconds.\n"), file=stderr()) } startpara[c("mu", "phi", "sigma")] <- as.list(apply(init.res$para, 2, median))
startlatent <- as.numeric(apply(init.res$latent, 2, median)) } if (length(mhcontrol) == 1) { cov.mh <- diag(mhcontrol, nrow = 4, ncol = 4) } else if (length(mhcontrol) == 2) { # in this case init.with.svsample > 0 # calculate proposal covariance matrix (in the order of svlsample) phi.t <- 0.5*log(2/(1-init.res$para[, "phi"])-1)
#rho.t we don't have
sigma2.t <- 2*log(init.res$para[, "sigma"]) mu.t <- init.res$para[, "mu"]
cov.mh <- cov(cbind(phi.t, rho.t = 0, sigma2.t, mu.t))
cov.mh[2, 2] <- mhcontrol$rho.var cov.mh <- mhcontrol$scale*cov.mh
} else {  # mhcontrol is a covariance matrix already
cov.mh <- mhcontrol
}

renameparam <- c("centered" = "C", "noncentered" = "NC")
if (!quiet) {
cat(paste("\nCalling ", asisprint(renameparam[parameterization], renameparam), " MCMC sampler with ", draws+burnin, " iter. Series length is ", length(y), ".\n",sep=""), file=stderr())
flush.console()
}

if (.Platform$OS.type != "unix") myquiet <- TRUE else myquiet <- quiet # Hack to prevent console flushing problems with Windows runtime <- system.time({ res <- svlsample_cpp(y, draws, burnin, designmatrix, thinpara, thinlatent, thintime, startpara, startlatent, priorphi[1], priorphi[2], priorrho[1], priorrho[2], 0.5, 0.5/priorsigma, priormu[1], priormu[2], priorbeta[1], priorbeta[2], !myquiet, myoffset, t(chol(cov.mh)), gammaprior, correct.latent.draws, parameterization, FALSE) }) if (any(is.na(res))) stop("Sampler returned NA. This is most likely due to bad input checks and shouldn't happen. Please report to package maintainer.") if (!quiet) { cat("Timing (elapsed): ", file=stderr()) cat(runtime["elapsed"], file=stderr()) cat(" seconds.\n", file=stderr()) cat(round((draws+burnin)/runtime[3]), "iterations per second.\n\n", file=stderr()) cat("Converting results to coda objects... ", file=stderr()) } colnames(res$para) <- c("mu", "phi", "sigma", "rho")
if (keeptime == "all") colnames(res$latent) <- paste0('h_', arorder + seq_along(y)) else if (keeptime == "last") colnames(res$latent) <- paste0('h_', arorder + length(y))
# create svldraws class
res$runtime <- runtime res$y <- y_orig
res$para <- coda::mcmc(res$para, burnin+thinpara, burnin+draws, thinpara)
res$latent <- coda::mcmc(res$latent, burnin+thinlatent, burnin+draws, thinlatent)
res$thinning <- list(para = thinpara, latent = thinlatent, time = keeptime) res$priors <- list(mu = priormu, phi = priorphi, sigma = priorsigma, rho = priorrho, gammaprior = gammaprior)
if (!any(is.na(designmatrix))) {
res$beta <- coda::mcmc(res$beta, burnin+thinpara, burnin+draws, thinpara)
colnames(res$beta) <- paste0("b_", 0:(NCOL(designmatrix)-1)) res$priors <- c(res$priors, "beta" = list(priorbeta), "designmatrix" = list(designmatrix)) } else { res$beta <- NULL
}
res$meanmodel <- meanmodel class(res) <- c("svldraws", "svdraws") if (!quiet) { cat("Done!\n", file=stderr()) cat("Summarizing posterior draws... ", file=stderr()) } res <- updatesummary(res, ...) if (!quiet) cat("Done!\n\n", file=stderr()) res } #' Minimal overhead version of \code{\link{svlsample}}. #' #' \code{svlsample2} is a minimal overhead version of \code{\link{svlsample}} #' with slightly different default arguments and a simplified return value #' structure. It is intended to be used mainly for one-step updates where speed #' is an issue, e.g., as a plug-in into other MCMC samplers. Note that #' absolutely no input checking is performed, thus this function is to be used #' with proper care! #' #' As opposed to the ordinary \code{\link{svlsample}}, the default values differ #' for \code{draws}, \code{burnin}, and \code{quiet}. Note that currently #' neither \code{expert} nor \code{\dots{}} arguments are provided. #' #' @param y numeric vector containing the data (usually log-returns), which #' must not contain zeros. Alternatively, \code{y} can be an \code{svsim} #' object. In this case, the returns will be extracted and a warning is thrown. #' @param 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 1. #' @param 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 0. #' @param priormu numeric vector of length 2, indicating mean and standard #' deviation for the Gaussian prior distribution of the parameter \code{mu}, #' the level of the log-volatility. The default value is \code{c(0, 100)}, #' which constitutes a practically uninformative prior for common exchange rate #' datasets, stock returns and the like. #' @param priorphi numeric vector of length 2, indicating the shape parameters #' for the Beta prior distribution of the transformed parameter #' \code{(phi + 1) / 2}, where \code{phi} denotes the persistence of the #' log-volatility. The default value is \code{c(5, 1.5)}, which constitutes a #' prior that puts some belief in a persistent log-volatility but also #' encompasses the region where \code{phi} is around 0. #' @param priorsigma single positive real number, which stands for the scaling #' of the transformed parameter \code{sigma^2}, where \code{sigma} denotes the #' volatility of log-volatility. More precisely, \code{sigma^2 ~ priorsigma * #' chisq(df = 1)}. The default value is \code{1}, which constitutes a #' reasonably vague prior for many common exchange rate datasets, stock returns #' and the like. #' @param priorrho numeric vector of length 2, indicating the shape parameters #' for the Beta prior distribution of the transformed parameter #' \code{(rho + 1) / 2}, where \code{rho} denotes the conditional correlation #' between observation and the increment of the #' log-volatility. The default value is \code{c(4, 4)}, which constitutes a #' slightly informative prior around 0 (the no leverage case) to boost convergence. #' @param thinpara single number greater or equal to 1, coercible to integer. #' Every \code{thinpara}th 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. #' @param thinlatent single number greater or equal to 1, coercible to integer. #' Every \code{thinlatent}th 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. #' @param thintime \emph{Deprecated.} Use 'keeptime' instead. #' @param keeptime Either 'all' (the default) or 'last'. Indicates which latent # volatility draws should be stored. #' @param quiet logical value indicating whether the progress bar and other #' informative output during sampling should be omitted. The default value is #' \code{TRUE}. #' @param startpara \emph{compulsory} named list, containing the starting values #' for the parameter draws. It must contain four #' elements named \code{mu}, \code{phi}, \code{sigma}, and \code{rho}, where \code{mu} is #' an arbitrary numerical value, \code{phi} is a real number between \code{-1} #' and \code{1}, \code{sigma} is a positive real number, and \code{rho} is #' a real number between \code{-1} and \code{1}. #' @param startlatent \emph{compulsory} vector of length \code{length(y)}, #' containing the starting values for the latent log-volatility draws. #' @return The value returned is a list object holding #' \item{para}{matrix of dimension \code{4 x draws} containing #' the \emph{parameter} draws from the posterior distribution.} #' \item{latent}{matrix of dimension \code{length(y) x draws} containing the #' \emph{latent instantaneous log-volatility} draws from the posterior #' distribution.} #' \item{meanmodel}{always equals \code{"none"}} #' @author Darjus Hosszejni \email{[email protected]@wu.ac.at} #' @seealso \code{\link{svlsample}} #' @keywords models ts #' @examples #' data(exrates) #' aud.price <- subset(exrates, #' as.Date("2010-01-01") <= date & date < as.Date("2011-01-01"), #' "AUD")[,1] #' draws <- svlsample2(logret(aud.price), #' draws = 10, burnin = 0, #' startpara = list(phi=0.95, mu=-10, sigma=0.2, rho=-0.1), #' startlatent = rep_len(-10, length(aud.price)-1)) #' @export svlsample2 <- function(y, draws = 1, burnin = 0, priormu = c(0, 100), priorphi = c(5, 1.5), priorsigma = 1, priorrho = c(4, 4), thinpara = 1, thinlatent = 1, thintime = NULL, keeptime = "all", quiet = TRUE, startpara, startlatent) { # Check whether 'thintime' was used if (!is.null(thintime)) stop("Argument 'thintime' is deprecated. Please use 'keeptime' instead.") if (keeptime == "all") thintime <- 1L else if (keeptime == "last") thintime <- length(y) else stop("Parameter 'keeptime' must be either 'all' or 'last'.") res <- svlsample_cpp(y, draws, burnin, matrix(NA), thinpara, thinlatent, thintime, startpara, startlatent, priorphi[1], priorphi[2], priorrho[1], priorrho[2], 0.5, 0.5/priorsigma, priormu[1], priormu[2], 0, 1, !quiet, 0, diag(0.1, nrow = 4, ncol = 4), TRUE, TRUE, rep(c("centered", "noncentered"), 5), FALSE) res$para <- t(res$para) res$latent <- t(res$latent) res$meanmodel <- "none"
rownames(res\$para) <- c("mu", "phi", "sigma", "rho")
res
}

.svsample <- function (...) {
.Defunct(new = "svsample2")
}

gregorkastner/stochvol documentation built on Sept. 3, 2019, 1:58 p.m.