article2.R
In stochvol: Efficient Bayesian Inference for Stochastic Volatility (SV) Models

## ----setup, include=FALSE, cache=FALSE------------------------------
knitr::render_sweave()
knitr::opts_chunk$set(prompt = TRUE,
                      fig.show = "hide",
                      warning = FALSE,
                      error = FALSE,
                      dev = c("pdf", "png"),
                      eval = FALSE,
                      dpi = 90,
                      message = FALSE,
                      echo = FALSE,
                      #cache = TRUE,
                      fig.path = "jss3813-",
                      tidy = FALSE)
base::options(continue = "+  ", prompt = "R> ", width = 70,
  useFancyQuotes = FALSE)

## ----presvrunmodel, eval=TRUE, results='hide'-----------------------
set.seed(1)
library("stochvol")
data("exrates")
ind <- which(exrates$date >= as.Date("2011-03-01") &
  exrates$date <= as.Date("2012-03-01"))
CHF_price <- exrates$CHF[ind]

## ----svrunmodel, results='hide', echo=TRUE--------------------------
#  set.seed(1)
#  library("stochvol")
#  data("exrates")
#  ind <- which(exrates$date >= as.Date("2011-03-01") &
#    exrates$date <= as.Date("2012-03-01"))
#  CHF_price <- exrates$CHF[ind]
#  res_sv <- svsample(CHF_price, designmatrix = "ar1")

## ----presvtrunmodel, results='hide', echo=TRUE, eval=TRUE-----------
set.seed(2)
CHF_logret <- 100 * logret(CHF_price)

## ----svtrunmodel, echo=TRUE, results='hide'-------------------------
#  set.seed(2)
#  CHF_logret <- 100 * logret(CHF_price)
#  res_svt <- svtsample(CHF_logret, designmatrix = "ar0")

## ----svlrunmodel, eval=TRUE, echo=TRUE------------------------------
set.seed(3)
X <- cbind(constant = 1, 100 * logret(exrates$USD[ind]),
  100 * logret(exrates$JPY[ind]))
res_svl <- svlsample(CHF_logret, designmatrix = X,
  thin = 20)

## ----svlplot, echo=TRUE, fig.height=4, results='hide', eval=TRUE----
plot(res_svl, showobs = FALSE, dates = exrates$date[ind[-1]])

## ----svlbetaplot, echo=2, fig.width=6.7, fig.height=3.5, eval=TRUE, results='hide'----
par(mar = c(2.5, 1.5, 0.5, 0.5), mfrow = c(3, 2), mgp = c(1.7, 0.5, 0))
for (i in seq_len(3)) {
  coda::traceplot(svbeta(res_svl)[, i])
  coda::densplot(svbeta(res_svl)[, i], show.obs = FALSE)
}

## ----printsummary, echo=TRUE, eval=TRUE, results="markup"-----------
summary(res_svl, showlatent = FALSE)

## ----svlpredict, echo=TRUE, eval=TRUE-------------------------------
set.seed(4)
pred_ind <- seq(tail(ind, 1), length.out = 5)
pred_X <- cbind(constant = 1, 100 * logret(exrates$USD[pred_ind]),
  100 * logret(exrates$JPY[pred_ind]))
pred_svl <- predict(res_svl, 4, newdata = pred_X)

## ----plotsvlpred, echo=2:4, eval=TRUE, fig.height=3.5, fig.width=10, results='hide'----
par(mgp = c(1.7, 0.5, 0))
obs_CHF <- 100 * logret(exrates$CHF[pred_ind])
qs <- t(apply(predy(pred_svl), 2, quantile, c(0.05, 0.5, 0.95)))
ts.plot(cbind(qs, obs_CHF), xlab = "Periods ahead", lty = c(rep(1, 3), 2),
  col = c("gray80", "black", "gray80", "red"))

## ----svroll, echo=TRUE, results='hide'------------------------------
#  set.seed(5)
#  res <- svsample_roll(CHF_logret, n_ahead = 1, forecast_length = 30,
#    refit_window = "moving", calculate_quantile = c(0.01, 0.05),
#    calculate_predictive_likelihood = TRUE)

## ----printpriordefault, echo=TRUE, results='hide'-------------------
#  svsample(CHF_logret, priormu = c(0, 100), priorphi = c(5, 1.5),
#    priorsigma = 1, priorbeta = c(0, 10000))
#  svtsample(CHF_logret, priormu = c(0, 100), priorphi = c(5, 1.5),
#    priorsigma = 1, priorbeta = c(0, 10000), priornu = 0.1)
#  svlsample(CHF_logret, priormu = c(0, 100), priorphi = c(5, 1.5),
#    priorsigma = 1, priorbeta = c(0, 10000), priorrho = c(4, 4))
#  svtlsample(CHF_logret, priormu = c(0, 100), priorphi = c(5, 1.5),
#    priorsigma = 1, priorbeta = c(0, 10000), priornu = 0.1,
#    priorrho = c(4, 4))

## ----printpriorspecdefault, echo=TRUE, results='hide'---------------
#  ps <- specify_priors(mu = sv_normal(mean = 0, sd = 100),
#    phi = sv_beta(shape1 = 5, shape2 = 1.5), rho = sv_constant(0),
#    sigma2 = sv_gamma(shape = 0.5, rate = 0.5), nu = sv_infinity(),
#    beta = sv_multinormal(mean = 0, sd = 10000, dim = 1),
#    latent0_variance = "stationary")
#  svsample(CHF_logret, priorspec = ps)

## ----misspecification, echo=TRUE, results='hide'--------------------
#  set.seed(3)
#  y <- svsim(50)$y
#  svsample(y, expert = list(correct_model_misspecification = TRUE))

## ----fsvprepdata, echo=TRUE, eval=FALSE, fig.width=10, fig.height=5, results='hide'----
#  library("factorstochvol")
#  library("zoo")
#  data("exrates", package = "stochvol")
#  m <- 6
#  n <- 200
#  y <- 100 * logret(tail(exrates[, seq_len(m)], n + 1))
#  y <- zoo(y, order.by = tail(exrates$date, n))
#  plot(y, main = "", xlab = "Time")

## ----preorder, echo=TRUE, eval=FALSE--------------------------------
#  preorder(y, factors = 2)

## ----findrestrict, echo=TRUE, eval=FALSE----------------------------
#  findrestrict(y, factors = 2)

## ----runmodel, echo=TRUE, eval=FALSE--------------------------------
#  set.seed(1)
#  res <- fsvsample(y, factors = 2, draws = 3000, zeromean = FALSE,
#    thin = 10, quiet = TRUE)

## ----printrres, echo = TRUE, eval=FALSE-----------------------------
#  res

## ----covn, echo = TRUE, eval=FALSE----------------------------------
#  dim(cov_n <- covmat(res))

## ----logdetcovn, echo = 2:5, fig.width = 10, fig.height=3.5, results='hide', eval=FALSE----
#  par(mfrow = c(1, 2), mgp = c(1.7, 0.5, 0), mar = c(3, 3, 1, 1))
#  logdet <- function (x) log(det(x))
#  logdet_n <- apply(cov_n[,,,1], 3, logdet)
#  ts.plot(logdet_n)
#  acf(logdet_n, main = "")

## ----covess, echo = TRUE--------------------------------------------
#  round(apply(cov_n, 1:2, coda::effectiveSize))

## ----corimageplot, echo=2:3, eval=FALSE, results='hide'-------------
#  par(mfrow = c(1, 3), xpd = TRUE)
#  corimageplot(res, these = seq(1, n, length.out = 3), plotCI = "circle",
#    plotdatedist = 2, date.cex = 1.1)

## ----voltimeplot, echo=2:3, eval=FALSE, fig.width = 10, fig.height = 3, results='hide'----
#  par(mgp = c(1.7, 0.5, 0), mar = c(2, 1.5, 1, 0.5))
#  palette(RColorBrewer::brewer.pal(7, "Dark2")[-5])
#  voltimeplot(res, legend = "top")

## ----cortimeplot, echo=2:4, eval=FALSE, fig.width = 10, fig.height = 5, results='hide'----
#  par(mfrow = c(2, 1), mgp = c(1.7, 0.5, 0), mar = c(2, 1.5, 1, 0.5))
#  palette(RColorBrewer::brewer.pal(6, "Dark2"))
#  cortimeplot(res, 1)
#  cortimeplot(res, 2)

## ----comtimeplot, echo=2, eval=FALSE, fig.height = 6.5, results='hide'----
#  par(mgp = c(1.7, 0.5, 0), mar = c(3, 3, 1, 1))
#  comtimeplot(res, maxrows = 6)

## ----loadplot2, echo=2:3, eval=FALSE, fig.width=4.5, fig.height=4.5, results='hide'----
#  par(mgp = c(1.7, 0.5, 0), mar = c(2.7, 2.7, 2, 0.5))
#  facloadpairplot(res)
#  facloadcredplot(res)

## ----varplot, echo=2, eval=FALSE, fig.width=10, fig.height=4, results='hide'----
#  par(mgp = c(1.7, 0.5, 0), mar = c(2.7, 2.7, 2, 0.5))
#  logvartimeplot(res, show = "fac")

## ----varplot2, echo=2, eval=FALSE, fig.width=7, fig.height=6.5, results='hide'----
#  par(mgp = c(1.7, 0.5, 0), mar = c(2.7, 2.7, 2, 0.5))
#  logvartimeplot(res, show = "idi", maxrows = 6)

## ----evdiag, fig.width=10, fig.height=4, eval=FALSE, echo=2:4, results = 'hide'----
#  par(mgp = c(1.7, 0.5, 0), mar = c(2.7, 2.7, 2, 0.5))
#  set.seed(6)
#  largemodel <- fsvsample(y, factors = 6)
#  evdiag(largemodel)

## ----predcov1, eval=FALSE, echo=TRUE--------------------------------
#  set.seed(4)
#  predcor1 <- predcor(res)
#  round(apply(predcor1[,,,1], 1:2, mean), 2)
#  round(apply(predcor1[,,,1], 1:2, sd), 2)

## ----preddist, eval=FALSE, fig.height = 6, fig.width = 9, echo = TRUE, results='hide'----
#  set.seed(5)
#  predcov_1 <- predcov(res)
#  preddraws <- t(res$beta)
#  for (i in seq_len(NROW(preddraws)))
#    preddraws[i,] <- preddraws[i,] + chol(predcov_1[,,i,1]) %*% rnorm(m)
#  plotlims <- quantile(preddraws, c(0.01, 0.99))
#  LSD::heatpairs(preddraws, labels = colnames(y), cor.cex = 1.5, gap = 0.3,
#    xlim = plotlims, ylim = plotlims)

## ----predloglik, eval=FALSE, echo = TRUE----------------------------
#  set.seed(6)
#  predloglik(res, matrix(0, nrow = 2, ncol = m), ahead = 1:2, each = 10)