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

#### Documented in svsim

#  #####################################################################################
#  R package stochvol by
#     Gregor Kastner Copyright (C) 2013-2020
#     Darjus Hosszejni Copyright (C) 2019-2020
#
#  This file is part of the R package stochvol: Efficient Bayesian
#  Inference for Stochastic Volatility Models.
#
#  The R package stochvol is free software: you can redistribute it
#  and/or modify it under the terms of the GNU General Public License
#  as published by the Free Software Foundation, either version 2 or
#  any later version of the License.
#
#  The R package stochvol is distributed in the hope that it will be
#  useful, but WITHOUT ANY WARRANTY; without even the implied warranty
#  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
#  General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with the R package stochvol. If that is not the case, please
#  #####################################################################################

#' Simulating a Stochastic Volatility Process
#'
#' \code{svsim} is used to produce realizations of a stochastic volatility (SV)
#' process.
#'
#' This function draws an initial log-volatility \code{h_0} from the stationary
#' distribution of the AR(1) process defined by \code{phi}, \code{sigma}, and \code{mu}.
#' Then the function jointly simulates the log-volatility series
#' \code{h_1,...,h_n} with the given AR(1) structure, and the log-return'' series
#' \code{y_1,...,y_n} with mean 0 and standard deviation \code{exp(h/2)}.
#' Additionally, for each index \code{i}, \code{y_i} can be set to have a conditionally heavy-tailed
#' residual (through \code{nu}) and/or to be correlated with \code{(h_{i+1}-h_i)}
#' (through \code{rho}, the so-called leverage effect, resulting in asymmetric log-returns'').
#'
#' @param len length of the simulated time series.
#' @param mu level of the latent log-volatility AR(1) process. The defaults
#' value is \code{-10}.
#' @param phi persistence of the latent log-volatility AR(1) process. The
#' default value is \code{0.98}.
#' @param sigma volatility of the latent log-volatility AR(1) process. The
#' default value is \code{0.2}.
#' @param nu degrees-of-freedom for the conditional innovations distribution.
#' The default value is \code{Inf}, corresponding to standard normal
#' conditional innovations.
#' @param rho correlation between the observation and the increment of the
#' log-volatility. The default value is \code{0}, corresponding to the basic
#' SV model with symmetric log-returns''.
#' @return The output is a list object of class \code{svsim} containing
#' \item{y}{a vector of length \code{len} containing the simulated data,
#' usually interpreted as log-returns''.}
#' \item{vol}{a vector of length
#' \code{len} containing the simulated instantaneous volatilities
#' \code{exp(h_t/2)}.}
#' \item{vol0}{The initial volatility \code{exp(h_0/2)},
#' drawn from the stationary distribution of the latent AR(1) process.}
#' \item{para}{a named list with five elements \code{mu}, \code{phi},
#' \code{sigma}, \code{nu}, and \code{rho}, containing
#' the corresponding arguments.}
#' @note The function generates the log-returns'' by
#' \code{y <- exp(-h/2)*rt(h, df=nu)}. That means that in the case of \code{nu < Inf}
#' the (conditional) volatility is \code{sqrt(nu/(nu-2))*exp(h/2)}, and that corrected value
#' is shown in the \code{print}, \code{summary} and \code{plot} methods.
#'
#' To display the output use \code{print}, \code{summary} and \code{plot}. The
#' \code{print} method simply prints the content of the object in a moderately
#' formatted manner. The \code{summary} method provides some summary statistics
#' (in \%), and the \code{plot} method plots the the simulated 'log-returns'
#' \code{y} along with the corresponding volatilities \code{vol}.
#' @author Gregor Kastner \email{gregor.kastner@@wu.ac.at}
#' @keywords datagen ts
#' @example inst/examples/svsim.R
#' @export
svsim <- function(len, mu = -10, phi = 0.98, sigma = 0.2, nu = Inf, rho = 0) {

name_len <- "length of simulated data set"
assert_numeric(len, name_len)
assert_single(len, name_len)
assert_ge(len, 2, name_len)

name_mu <- "input parameter mu"
assert_numeric(mu, name_mu)
assert_single(mu, name_mu)

name_phi <- "input parameter phi"
assert_numeric(phi, name_phi)
assert_single(phi, name_phi)
assert_gt(phi, -1, name_phi)
assert_lt(phi, 1, name_phi)

name_sigma <- "input parameter sigma"
assert_numeric(sigma, name_sigma)
assert_single(sigma, name_sigma)
assert_positive(sigma, name_sigma)

name_nu <- "input parameter nu"
assert_numeric(nu, name_nu)
assert_single(nu, name_nu)
assert_gt(nu, 2, name_nu)

name_rho <- "input parameter rho"
assert_numeric(rho, name_rho)
assert_single(rho, name_rho)
assert_gt(rho, -1, name_rho)
assert_lt(rho, 1, name_rho)

len <- as.integer(len)

h <- rep_len(as.numeric(NA), length.out=len)
h0 <- rnorm(1, mean=mu, sd=sigma/sqrt(1-phi^2))
tau <- if (is.finite(nu)) {
1/rgamma(len, shape=nu/2, rate=nu/2-1)
} else {
rep_len(1, length.out=len)
}
eta <- rnorm(len)
eps <- rho*eta + sqrt(1-rho^2)*rnorm(len)

# simulate w/ simple loop
h <- mu + phi*(h0-mu) + sigma*rnorm(1)  # same marginal distribution as h0
for (i in seq_len(len-1)) {
h[i+1] <- mu + phi*(h[i]-mu) + sigma*eta[i]
}
y <- exp(h/2) * sqrt(tau) * eps  # "log-returns"

ret <- list(y = y,
vol0 = exp(h0/2),
vol = exp(h/2) * sqrt(tau),
para = list(mu = mu,
phi = phi,
sigma = sigma,
rho = rho,
nu = nu),
latent = h,
latent0 = h0,
tau = tau)
class(ret) <- "svsim"
ret
}

#' @export
print.svsim <- function(x, ...) {
cat("\nSimulated time series consisting of ", length(x$y), " observations.\n\n", "Parameters: level of latent variable mu = ", x$para$mu, "\n", " persistence of latent variable phi = ", x$para$phi, "\n", " standard deviation of latent variable sigma = ", x$para$sigma, "\n", " degrees of freedom parameter nu =", x$para$nu, "\n", " leverage effect parameter rho =", x$para$rho, "\n", sep="") cat("\nSimulated initial volatility:", x$vol0*x$correction, "\n") cat("\nSimulated volatilities:\n") print(x$vol*x$correction, ...) cat("\nSimulated data (usually interpreted as 'log-returns'):\n") print(x$y, ...)
}

#' @export
plot.svsim <- function(x, mar = c(3, 2, 2, 1), mgp = c(1.8, .6, 0), ...) {
op <- par(mfrow = c(2, 1), mar = mar, mgp = mgp)
plot.ts(100*x$y, ylab = "", ...) mtext("Simulated data: 'log-returns' (in %)", cex = 1.2, line = .4, font = 2) plot.ts(100*x$vol, ylab = "", ...)
mtext("Simulated volatilities (in %)", cex = 1.2, line = .4, font = 2)
par(op)
}

#' @export
summary.svsim <- function(object, ...) {
ret <- vector("list")
class(ret) <- "summary.svsim"
ret$len <- length(object$y)
ret$para <- object$para
ret$vol <- summary(100*object$vol)
ret$y <- summary(100*object$y)
ret$vol0 <- 100*object$vol0
ret
}

#' @method print summary.svsim
#' @export
print.summary.svsim  <- function(x, ...) {
cat("\nSimulated time series consisting of ", x$len, " observations.\n", "\nParameters: level of latent variable mu = ", x$para$mu, "\n persistence of latent variable phi = ", x$para$phi, "\n standard deviation of latent variable sigma = ", x$para$sigma, "\n degrees of freedom parameter nu = ", x$para$nu, "\n leverage effect parameter rho = ", x$para$rho, "\n", sep="") cat("\nSimulated initial volatility (in %): ") cat(x$vol0, "\n")
cat("\nSummary of simulated volatilities (in %):\n")
print(x$vol) cat("\nSummary of simulated data (in %):\n") print(x$y)
invisible(x)
}


## Try the stochvol package in your browser

Any scripts or data that you put into this service are public.

stochvol documentation built on July 12, 2021, 5:08 p.m.